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VARIANTS AND EXTENSIONS OF A FAST DIRECT NUMERICAL CAUCHY-RIEMANN 


SOLVER, WITH ILLUSTRATIVE APPLICATIONS 
E. Dale Martin and Harvard Lomax 
Ames Research Center 


SUMMARY 


Revised and extended versions of a fast, direct (noniterative) numer- 
ical Cauchy-Riemann solver are presented for solving finite-difference 
approximations of first-order systems of partial differential equations. 
Although the difference operators treated are linear and elliptic, one signif- 
icant application of these extended direct Cauchy-Riemann solvers is in the 
fast, semidirect (iterative) solution of fluid-dynamic problems governed by 
the nonlinear mixed elliptic-hyperbolic equations of transonic flow. Different 
versions of the algorithms are derived and the corresponding FORTRAN computer 
programs for a simple example problem are described and listed. The algorithms 
are demonstrated to be efficient and accurate. 


INTRODUCTION 


Highly efficient finite-difference algorithms and programs called "fast 
direct elliptic solvers" are becoming increasingly more useful for solving 
partial differential equations in practical problems as the techniques that 
incorporate the solvers, as well as the solvers themselves, become more highly 
developed. The fast direct solvers differ from other elliptic solution algo- 
rithms in that they take advantage of the sparseness and regularity of the 
highly ordered coefficient matrix of the set of finite-difference equations 
to obtain a direct solution in a sequence of simple recursive operations. 

Thus, they are essentially noniterative, and the consequent high efficiency 
accounts for their attractiveness in various applications. 

The purposes of this report are (a) to develop highly efficient and 
stable numerical algorithms called Cauchy-Riemann solvers for use in methods 
that treat the first-order elliptic operators in generalized Cauchy-Riemann 
equations, (b) to describe the development of several different versions of 
these algorithms, and (c) to provide a simple example computer program for 
each version, along with the subroutines for the Cauchy-Riemann solvers. 

Cauchy-Riemann solvers have significant applications in semidirect 
(globally implicit) iteration techniques for rapid solution of transonic 
flows. It is expected that, as the direct solvers become further developed 
and generalized, the significance of their applications will increase. 



Fast direct elliptic solvers were first developed for solving Poisson T s 
equation on a rectangle (refs. 1-5). Since then, the algorithms have been 
generalized, extended, and applied in numerous ways: (a) Fast solvers have 

been extended to three dimensions (e.g., refs. 5 and 6, among others). 

(b) They have been further developed and extended to treat more general 
second-order elliptic equations and more general boundary conditions and 
meshes (refs. 5 and 7-17). (c) They have been extended to include interior 

conditions and irregular domains (refs. 18-23). (d) They have been extended 

to fourth-order (biharmonic) partial differential equations in references 22 
and 24. (e) They have been extended to apply to sets of first-order elliptic 

equations in reference 25. Chapter 9 in reference 26 also describes an early 
direct solver for the Cauchy-Riemann equations. However, reference 26 indi- 
cates that the algorithm described there has a very low mesh-number limitation, 
which would severely restrict the usefulness; in addition, that algorithm 
would be significantly less efficient than the "fast" solvers of references 1 
to 5 and 25 that are based on either fast Fourier transforms or cyclic reduc- 
tion. (f) Direct elliptic solvers have been applied within iteration schemes 
(i.e., semidirect methods) to extend their usefulness to nonseparable elliptic 
equations (refs. 27 and 28) and either to Poisson equations as part of a 
system of nonlinear equations (refs. 29 and 30; see also ref. 31) or as the 
total driving algorithm in a semidirect iteration of nonlinear equations 
(refs. 32-34) and of equations that are both nonlinear and nonelliptic in 
some regions of the solution domain (refs. 35-38). (g) In addition to other 

applications in the literature (now too numerous to attempt to list here) , a 
notable application of a direct Poisson solver is as an iterative-acceleration 
device for finite-difference, line-relaxation solutions of the nonlinear, 
mixed-type equations of transonic flow (refs. 39-41). 

The fast Cauchy-Riemann solver developed in reference 25 and used in 
reference 35, and the more recent variant and extension of it used in refer- 
ences 36 to 38 (to be described below as version C), can be used in problems 
where a Poisson solver or other second-order elliptic solver could also be 
used. In addition, however, the direct solution of the corresponding set of 
first-order equations can be especially useful, for example, in fluid dynamics, 
because (a) the boundary conditions, including solid-surface conditions and 
and internal jump (e.g., shock-wave discontinuity or contact discontinuity) 
conditions, are given most naturally in terms of components of the velocity 
vector; and (b) the formulation in terms of velocity components (generalized 
Cauchy-Riemann equations) allows both point sources and point vortices to be 
included simultaneously. This is contrasted with stream-function and velocity- 
potential formulations, which are second-order partial differential equations 
(where velocity components are defined in terms of derivatives of either 
stream function or velocity potential). The stream-function formulation allows 
specification of point vortices but not point sources of mass, and the 
velocity-potential formulation allows specification of point sources but not 
point vortices. 

In the sections below, after description of the computational meshes and 
mathematical definitions, five versions of the Cauchy-Riemann solver, denoted 
A to E, are developed. In addition, computer programs including all subrou- 
tines for each version are listed in appendices for a simple example problem 
that illustrates use of the solvers. Version A is the original version of 
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the Cauchy-Riemann solver derived in reference 25 and used for the slightly 
supercritical transonic flow problem in reference 35. The present matrix 
derivation of version A is more concise than the previous derivation, and 
more details including the program listing can be given here. A significant 
step for improving the efficiency resulted in the development of version B. 
Versions B and C are essentially the same except that additional terms (with 
constant coefficients) needed to stabilize iterations in a semidirect method 
are included in C. Version C is the solver used for the transonic-flow cal- 
culations in references 36 to 38. Versions D and E are developed for future 
use. They use a mesh that is symmetrical about the x axis, whereas the 
previous versions used different configurations for upper and lower boundaries 
on the mesh and different conditions applied there. These versions are 
especially desired for future treatment of airfoils represented by interior 
conditions specified on a slit in the center of the mesh. Version D uses 
constant coefficients on the extra stabilizing terms, whereas version E uses 
arbitrarily specified, variable (in x) coefficients on the extra terms. 

Although versions B and C together and D and E together are derived as single 
algorithms, of which the pairs of versions are special cases, the separate 
sample programs are provided for all five versions because of the differences 
in efficiency that may be significant in future applications. The distinguish- 
ing features of each version are further described in the sections dealing 
with the derivations (beginning on pages 18, 21, and 24). A guide for the 
prospective user’s choice of algorithm is given at the end of the section 
preceding the Concluding Remarks. 


GENERAL DIFFERENTIAL AND FINITE-DIFFERENCE EQUATIONS 


The algorithm and programs to be developed are for numerical solution of 
the nonhomogeneous Cauchy-Riemann equations: 


dv 


s(x 9 y) 


(la) 


dU _ %V_ 

dx 


-u(x 9 y) 


(lb) 


or of more general extensions of equations (1), which have quite general 
physical applications. Throughout this report, the terminology of fluid 
mechanics is used. The dependent variables u and V then frequently repre- 
sent velocity components. Generally, however, for any two-dimensional vector 1 
(u,y), equation (la) is the definition of the divergence (a) of the vector, 
and equation (lb) is the definition of the curl (co) of the vector or, equiva- 
lently, the "vorticity" of the vector (u 9 v). Any point where the divergence 
s is not zero is called a "source" of the vector (u 9 v) , and any point where 
a) is not zero is called a "vortex" of the vector (u 9 v). In any region where 


1 A11 equations and definitions in this section have three-dimensional 
counterparts, but the present scope is limited to two dimensions. 
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s is zero, the vector ( u,v ) is "solenoidal and in any region where to is 
zero, the vector ( u,V ) is "irrotational . " Consistent with the terminology of 
fluid mechanics and the above general definitions, then, s(x,y) is here called 
a "source function," to (x,y) is called a "vorticity function," and equations (1) 
are referred to, respectively, as the "continuity equation" and the "rotation- 
ality equation." 

The algorithms to be developed can be useful in the efficient numerical 
solution of generalized Cauchy-Riemann equations; that is, this work is 
relevant in certain ways to the following more general form of the set of 
differential equations : 


lx + dy ~ s ( x >y> u > v ’* x > d y ') (2a) 

Hy ~ ~dx ~ - t *( X ’y> U > V ’ d x ’ d y') (2b) 


where the notation on the right side indicates that s and oj can be arbi- 
trarily specified functions not only of x and y but also of the dependent 
variables, u and V , and derivatives of arbitrary order of u and V with 
respect to x and y , with being the partial-derivative operator, etc. 

Thus the equations can be both nonlinear and nonelliptic, even though the left 
side of equations (2) is a linear elliptic operator. 


Iterative schemes for the solution of equations (2) treat the right side 
as known from a previous iteration. At each iteration, then, the "iteration 
equations" for solving the system of equations (2) are actually of the form 
of equations (1), in the simplest case, or they may be of a slightly more 
general linear form. 


The algorithms to be derived for solving equations (1) use central finite 
differences on the left side. For use in iteration schemes, it has been found 
(refs. 36-38) that, in some cases, certain stabilizing terms need to be added 
to both sides of the equations. Then the difference operator on the left will 
contain terms in addition to those representing the derivatives in equa- 
tions (1) . The most general form of the approximate set of finite-difference 
equations to be considered here for representing equations (1) is 


u 


+ 


- u 


X 


+ 


- 

h 


+ a + u + + a u 


* + 
- vT 

y 



a + u + + a u + s ( • ) 
o o 

s' (•) 


- 0 )(+) 


(3a) 

(3a T ) 
(3b) 


These equations are written for a staggered mesh, which has been found to have 
numerous advantages. Refer to figure 1(a), on which: the point at which 

equation (3a) is being solved is represented by (•); the point at which equa- 
tion (3b) is being solved is represented by (+) ; the values of u and V at 
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the centers of the boundaries of the mesh cell for the continuity equation 
(centered at the dot) are u + , u“, V + 9 and V ~ 9 and the values of u and V 
at the centers of the boundaries of the mesh cell for the rotationality 
equation (centered at the cross) are u * , , V* , and The notations s ( • ) 

and co (+) in equations (3) indicate that the source function s is evaluated 
at the dot and the vorticity function co is evaluated at the cross. The 
quantities a + and a" in equation (3a) may generally vary with x. The 
symbols uj and Uq in equation (3a) denote some known representations of 
ifr and u“, such as an analytical solution or results from a previous itera- 
tion in an iterative solution. The quantity s’(*) is defined to equal the 
right side of equation (3a). The x and y mesh intervals, h x and hy , are 
constants . 


COMPUTATIONAL MESHES AND DISCRETIZED VARIABLES 


For the development of the computational algorithms to solve equations (3), 
the configuration of mesh cells shown in figure 1(a) is imbedded in a large 
computational mesh. The different versions of the Cauchy-Riemann solver to 
be described use different treatments of the upper boundary, and so the mesh 
configurations differ there. However, since the left and bottom boundaries 
are the same for all versions, a system of mesh-point indexing can be used 
that is common to all versions, with origins of coordinate indices at or near 
the left bottom corner. 


h x 




T 

hy 

1 


(b) Indexing of basic configuration. 


Figure 1.- Mesh cells for continuity and rotationality equations. 
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Consider rectangular coordinates x and y and let the left and right 
boundaries be x = Xi and x = x u , and the bottom and top boundaries be y - y% 
and y = y u . Let 3 and k index the x and y directions, respectively, so 
that discrete coordinates Xj and yy_ are given by 

X q = x l + yh x * y k = y i + 

where J and k are integers. Now consider the indexing of the variables in 
figure 1(a) and refer to figure 1(b). For integers j and k 9 the point where 
oo is defined for equation (3b) is labeled oo^ ^ on figure 1(b). However, 
the point where the source function s is defined is displaced from 9 k 
by half intervals. To avoid the use of indices containing fractions (half 
intervals), define the discrete coordinates as 

x j' = x i + ( c ' ~ ~i) h x ’ y k ' = y i + (*' ~ 2 ) h y (4b) 

where j f and & f are integers. In order that the indices on oo have the 
same values as those on s corresponding to the rotationality equation and 
the continuity equation, respectively, for the staggered mesh cells in 
figures 1, let 

r = 3 , *' = k + \ (5) 


The locations where u and v are defined in figure 1(a) make it convenient 
to then define 

u i,k '’ u+ • v r,k’ v+ <6) 

We thus define values of to, s, u, and V at discrete points in figure 1(b) as 


*j,k “ (a? j ,y k) * 

k 1 = u ’ y k' ^ 9 


V j,k = v(x }^ y k ) ) 


(7) 


As used in equation (3a), the second-order-accurate, central finite-difference 
approximations needed for the continuity equation are then 


( du/ 3x) ,y * f 
d » ^ 


(u j,r u j-i,k' )/h x 


(3u/3 y)j, tk , 


{V j\k V o\k-i )/h y 


(8a) 
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Figure 2 Staggered computational meshes 








and, as used in equation (3b), the second-order-accurate, central finite- 
difference approximations needed for the rotationality equation are 

} (8b) 

* <*>+* ,k - v i\k )/h x ) 

The motivation for using integer values of j, j T , k 9 and k' in equations (4) 
and (7) is that, for purposes of indexing the variables u 9 V 9 s 9 and to, there 
need not be any distinction between j and J 1 or between k and k ' . Thus, 
for example, whenever uj ^ is written, uj ^ is implied. The distinction 
is important only when determining actual locations on the mesh. 

Figure 2 shows three different mesh configurations that correspond to 
the algorithm versions indicated. Figure 1(b) is imbedded in each mesh, as 
indicated by crosshatching (with j ,k = 2,2). The three meshes differ only 
at the upper boundary; the reasons for the differences are explained in later 
sections. The circled symbols in figure 2(a,b,c) indicate quantities that are 
specified as boundary conditions. As in figure 1(a), the dots in figure 2 
indicate points where the continuity equation (3a) is to be satisfied (values 
of s j\k' specified) and the crosses indicate points where the rotationality 
equation (3b) is to be satisfied (values of 9 k specified). An exception 
to this is that on the line y - y u in figure 2(a) for version A the 
n circled u" symbols are at points that should also be marked by crosses 
because the rotationality equation is to be satisfied there and values of a) , 
as well as u 9 are to be specified there. 

For defining the mesh dimensions and the dimensions of relevant matrices, 
let m 9 77? i, n 9 n\ 9 and n 2 be integers so that: 

m = the number of discrete x values where either the continuity or the 

rotationality equation is to be satisfied, that is, the number of dots or 
crosses in a horizontal row in figures 2(a), (b) , and (c) (all versions 
of the solver); 

has the value 4 for the special case illustrated in figure 2; 
mi denoted as the "mesh number for the x direction”: 

= m + 1; 

has the value 5 for the special case illustrated in figure 2; 

n = number of discrete y values where the rotationality equation is to be 
satisfied in figures 2(b) and 2(c), that is, the number of crosses in a 
vertical column in figures 2(b) and 2(c) (versions B, C, D, and E of the 
solver) ; 
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= also the number of discrete y values where the continuity equation is 

to be satisfied in figure 2(b), that is, the number of dots in a vertical 

column in figure 2(b) (versions B and C of the solver); 

has the value 3 for the special case illustrated in figure 2; 

n\ denoted as the "mesh number for the y direction"; 

= n + 1; 

= 2^, where L is a positive integer; 

= number of discrete y values where the rotationality equation is to be 
satisfied in figure 2(a), that is, the number of crosses in a vertical 
column in figure 2(a) (version A); 

= also the number of discrete y values where the continuity equation is 
to be satisfied in figures 2(a) and 2(c), that is, the number of dots in 
a vertical column in figures 2(a) and 2(c) (versions A, D, and E) ; 

has the value 4 for the special case illustrated in figure 2; 

Yl2 ~ vi + 2 ; 

has the value 5 for the special case in figure 2(a). 

To summarize the ranges of the discretized functions on the meshes in 
figure 2, table I is provided. 


TABLE I . - RANGES OF DISCRETIZED FUNCTIONS 





Values of coordi- 

Discrete 

Fun r f 1 on 

Version 

Ranges of 

nate where func- 

coordinate 


of solver 

indices 

tion is defined 

equation 

U 3>k' 

All 

j - 0 to m 

x = x . 

X j = X l + Jh x 


All 

k' = 1 tO Yi\ 

y = y' v 

y k ' = y SL + (k ' " h)h y 


A 

k' = n 2 

y = y u 

y u =y l + n '\ 

V j',k 

All 

j’ = 1 to mi 

X = 

X j' = x i + ^ ' ~ h)h x 


A, L>,E 

k - 0 to n\ 

y = y k 

y k = y i + kh y 


B,C 

k ~ 0 to n 

y = y k 

y k = y l + kh y 

s rx 

All 

j' - 1 to m 

X = 

x r = x n + (J ” ' h)h x 


A,D,E 

k 1 - 1 to Yl\ 

y = 

y r = vt + 


B,C 

k' = 1 to n 

y = y k’ 

V'v = + (fe ’ - h)h y 


All 

j - 1 to m 

x = x . 
3 

x , = x n + jh 
3 l x 


A 

k - 1 to ni 

y = y k 

y, = y n + kh 
y k v l y 


B,C,D,E 

k = 1 to n 

y = y k 

w T = y n + kh 
^k v l y 
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In terms of the above defined parameters, the mesh intervals h x and hy 

are 

h - (x - x )/[ m + (1/2)] , all versions (9a) 

CC U Jo 

h = (y u ~ y^/n 1 , versions A,D,E (9b) 

h = (U u ~ y^ln + (1/2) ] , versions B,C (9c) 

It will be convenient for the FORTRAN programs to also define "nominal" values 
of x u , x^, and y u , denoted as x un , x ln , y un , by 


h = (x - x. ) /m-> 
x un In 1 


\ ■ Kn - V /n l 


(9d) 

(9e) 


The usefulness of these definitions is explained in the later section "FORTRAN 
Programs . " 


MATRIX DEFINITIONS AND RELATIONS 


With the above definitions of mesh parameters and variables, we can now 
define some particular matrices and consider some matrix relations that will 
be of common use in the next four sections for developing the different ver- 
sions of the Cauchy-Riemann solvers. As described above, we do not distin- 
guish between J and j T or between k and k 1 for purposes of indexing the 
dependent variables and specified functions. 


First define the column vectors of dimension m: 


C ° 1[ “l ,k’ U 2,k ’ ■ 

. . , u 7 ] 
’ m,k 

(10a) 

c ° 1[t \ ,v v 2,r ■ 

. . , v 7 ] 
’ m,k 

(10b) 

o 

o 

H 

N> 

\0 

* * ^m,k) 

(10c) 

col[ »i ,*• »2 ,r ■ 

* •’ g m,k ] 

(10d) 


where fj y k anc ^ 9 k are to ke defined later in terms of Sj 9 k an d >k* 
The range of k in equations (10) is slightly different for different ver- 
sions of the solver and is specified later. 
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With the vectors in equations (10) as elements, next define the following 
block column vectors of block dimension p or q: 


U = col[U l5 U 2 , . . . , U ] 


(Ha) 


V = col[Vi, V 2 V ] 

r 


(lib) 


f = col[f ! , f 2 , . . f ] 

H 


(lie) 


g = coltg!, g 2 , . . . , g^] 


(lid) 


In each application of these vectors, p and q will be specified in terms of 
n and n\. 

A general tridiagonal matrix of square dimension m is denoted by 

r*i i 


T (a.,b . 9 q .) = 
m 0 3 Q 


a . b . c . 

0 J 3 


a b 
m m 


Denote a general diagonal matrix of square dimension m by 


D (b .) = T (0 9 b .,0) 
m 3 m 9 j* 


Special cases are the unit matrix and null matrix, both of dimension m: 


I = 1^(0, 1,0) , O = 1^(0, 0,0) 


Equations (12) and (13) may also be used to define general block-tridiagonal 
and block-diagonal matrices of block dimension p, for example: 


I 



T (A,B,C) = 


B 

A 


C 

B 


ABC 
A B 


(15a) 


pxp blocks 


and 


D (B) = T (0,B,0) 

p p 


(15b) 


where the arguments A, B, and C are dummy m*m matrices and 0 is defined in 
equation (14). At this point, it is also convenient to define the following 
special partitioned square block matrices of block dimension p: 


H = 

P 


J (B) = 
P 


"o 


0 

1 

1 

I ’ 

0 

• • • 

0 

1 

I 

I 

• 


• 

1 

• 

• 


• 

1 

• 

_0 

• . . 

0 

1 

1 

1 

I _ 

"o 

... 

0 

1 

o" 



: 

1 

• 

• 


• 

1 

• 

0 


0 

1 

1 

0 

_0 


0 

1 

1 

B _ 


-* pxp blocks 


' pxp blocks 

and in terms of equations (15b) and (17) then define 


(16) 


(17) 


D'(B) - D p (B) - j J p (B) 


(18) 


which is a block-diagonal matrix with the last diagonal block equal to half the 
matrix occupying each of the other diagonal blocks. 

A convenient matrix for following developments is denoted by Ep^(A) and 
is defined as follows. Consider a general matrix M of dimension p x a: 

pq * H 


M 


pq 


mn 

™ 12 

^13 

... 

^1 q 

m 2 1 

m 22 




: 

■ 

. 



• 

• 




m p\ 

m p2 



1 

Cr 


(19) 
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Let l = min {p 9 q}. Then the diagonal is defined to consist of the elements 
m ii> i = 1, 2, . . . , £ (whether or not p = q; see, e.g., ref. 42); the 
superdiagonal consists of the elements m j and the subdiagonal consists 

of the elements m i 9 i-\* Now to define Ep^(A), let each element on its 
diagonal be a matrix -A of dimension m and let each element on its super- 
diagonal be the matrix A. Let all other elements be ?77-dimensional null 
matrices. Thus Ep^(A) has p rows and q columns of ^-dimensional blocks. 
Two special cases of interest are for q = p + 1: 


E 


P>P + 1 


(A) 


-A A 

O -A A 


and for q = p: 


O 


-A 


A 


px(p+ l) blocks 


-A A 


E (A) E E (A) 

P VP 


0 -A 


A 


O 


-A 


J pxp blocks 


(20a) 


(20b) 


It is convenient to further define some special matrices that will be of 
use in the following four sections. From here on, whenever one of these 
matrix symbols is used, it has the specific definition given in this section. 
For these definitions, let 3, an d 32, j be quantities that are to be 

specified later. Then, in terms of the above-defined notation, define the 
particular w-dimensional square matrices: 


and 


AS 


( 21:0 


B = T m (o,e,-e) 


(21b) 


C = 21 + AB 

= T (-a.,B.,-y.) (21c) 

m o 3 3 
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where 




a . 
3 


33 0 A 


6 i ■ 2 + B<B 2,i + hj 


(j = 1 to m) 


yjt ■ 


BB .,i 


(21d) 


with 82 1 = 0 ; also define the particular block matrices of block dimension 
p, or q 9 or p x q 9 as indicated: 



$ 

II! 

(22a) 

d b 

03 

til 

(22b) 

d 'b 

= D' (B) 

P 

(22c) 

E 

III 

(22d) 

T 

E 1 

= Block transpose of E 

(22e) 


where A and B are defined in equations (21) and p and q are to be specified 
in each application. 

The following relationships are then easily found and will be of particu- 
lar use, for either q = p or q = jp + 1: 



M 

II 

> 

(23a) 


[D (A) ]E = E (A) 

p Ph 

(23b) 


[Pp (A) ] [D (B) ] = D p (AB) 

(23c) 


[D p (A)][DpB)] = D p (AB) 

(23d) 

For q = p: 

EE r = T (- 1 , 21 , -I) - J (I) 

(23e) 

but, for q = p + 1 

: 



EE T = T p (-I,2I,-I) 

(23f) 

As a consequence of 

equations (23a) and (23b) , we find for either 

q = p or 


q = p + 1: 
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I 


ED a = [D (A)]E (24a) 

and from equations (23c) to (23f) , along with (21c), that for q = p: 

EE T + [D p (A)]D fl = T p (-1, 21,-1) - J (I) + D p (AB) 

= T p (-I,C,-I) - J (I) (24b) 

EE T + [D p (A)]D^ = T p (— 1,21,— I) - J p (I) + D^(AB) 

= T p (-I,0,-I) + D^(C) (24c) 

and for q = p + 1: 

EE T + [D (A)]D b = T (-I,C,-I) (24d) 

SIMPLEST FORMULATION OF MATRIX EQUATIONS, AND MOTIVATION 
FOR ALTERNATIVE VERSIONS 


Consider as the simplest and most straightforward formulation of equa- 
tions (3) the case where the mesh is as shown in figure 2(b) except that the 
circles are removed from the row of circled V symbols in the mesh row 
k ~ n. That is, values of u are specified at k' - (y = y u ) , but v is 
not specified at k = n. 

Define 


3 = h /h (25) 

y x 

and assume for the formulation in this section that a + and oT in equa- 
tions (3) are zero. Then, with the definitions in equations (4) to (8), 
equations (3) become 


e(u j,k 


- u 




V j,k-1 \ S j,k 


-u 


i,k + u j,k + i + B(u 


- v . , , 7 ) - -h co . 7 

<7 ,k j+i ,k y j 9 k 


(26a) 


(26b) 


for all j - 1 to m and k = 1 to n. (The value of m is arbitrary, but 
ri\ = n + 1 - 2^ where L is an integer. This is important for the discussion 
of the most efficient cyclic reduction of the equations that result later.) 
Equations (26) are second-order accurate. With all specified boundary values 
of u j 9 k anc * v j 9 k ( as indicated in fig. 2(b), except as noted above) moved to 
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the right sides of equations (26), denote the resulting right side of equa- 
tion (26a) by fj 9 k anc * t * ie resulting right side of equation (26b) by 
Then with the vector and matrix definitions in the above section and with 


and 


= = 3 (except $ 2>1 = 0) 


( 27 ) 


p = q = n (28) 

the complete set of finite-difference equations (26) for all j 9 k can be 
written as 

n blocks n blocks 


(29) 


or, more concisely, 

D^U - E 7 V = f (30a) 

EU + D^V = g (30b) 

Thus, equation (30a) represents the entire set of continuity equations (26a), 
and equation (30b) represents the entire set of rotationality equations (26b). 

Equations (30) can be reduced by eliminating U. Premultiply equa- 
tion (30a) by -E, premultiply equation (30b) by D n (A) , add the two resulting 
matrix equations, and use the relations (24a) and (24b) to obtain 


~A 




1 

I 

0 




u i 




A 



1 

1 

1 

-I 

I 

• 

• 



u 2 


*2 





1 




o 


# 







1 





* 


• 




A 

1 



-I 

I 


u n 



- - 

- - - 

- - - 

- - 

4- 

- - 

- - - 

- - - 

- - - 


- “ 

= 

- ~ 

-I 

I 



1 

1 

B 





Vi 


§1 

0 

-I 

x 


1 

1 


B 




V 2 




# . 

x 

I 

1 

1 



*. 



• 

• 


• 



0 

-I 

1 

1 




B 


JV 


_ §n_ 



where 


F = col[F 1} F 2 , . . *» F p] (32a) 

= [D (A)]g - Ef (32b) 

or 

f k • As k + h - '*+! > k * 1 to p <32c) 

where, in this case, p = n, and f must be defined to be zero when p = q 

for equation (32c) to apply; and where the Q component of the vector repre- 
sented by each product A is just (j s 1 to w, 

k = 1 to p) in which 3 2 ,1 = 0 as in equations (27). 

Equation (31a) is a block- tridiagonal equation to solve for V, with the 
right side, F, easily obtained from the simple operations indicated in equa- 
tion (32c). In an expanded form, the block-tridiagonal equation (31a) is 


"C -I 
-I 

c 

-I 

Equation (31b) could be solved (ref. 25, p. 68) by a fast, direct elliptic 
solver for each v j >k> and then the bottom half of equation (29) could be used 
to solve for each However, the direct solution of equation (31b) is 

significantly more difficult than it would be if the entry in the nth column 
of the nth row of the block coefficient matrix were C, like all the other 
diagonal entries , rather than C - I . 


-1 

C-I 



-vr 


"Pi" 


V 2 


F 2 


• 


• 


• 


• 


V 


F 


l n -* 


L 71 J 


(31b) 


There are two approaches that one can now take. Equation (31b) is in a 
form equivalent to that obtained for a Poisson equation with a Neumann condi- 
tion specified in a certain way at one boundary (which would result in the 
last diagonal entry, C-I). Therefore, one possible approach is that used by 
Sweet (ref. 9) and by Schumann and Sweet (ref. 17) in solving Poisson’s equa- 
tion: to deal directly with the irregular element of the block coefficient 
matrix by developing a separate and more complicated factorization for calcu- 
lations during the reduction process involving that element. A different point 
of view avoids the additional complication both in the algorithm development 
and in the programming. This approach, which had been taken in reference 25 
and also is taken here, is to look for simple ways of modifying the problem 
formulation so as to obtain a very regular coefficient matrix, with all block- 
diagonal elements the same, so that simpler algorithms and programs can be 
used. The algorithm versions described in the next three sections deal with 
the alternative formulations. 
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VERSION A - ORIGINAL FAST DIRECT CAUCHY-RI EMANN SOLVER 


Version A described here is the original version of the Cauchy-Riemann 
solver developed in reference 25 and used for fast iterative calculations in 
a nonlinear problem in reference 35. The derivation is given here, not only 
for completeness in the sequence of development both of algorithms and of 
programs, but also because the matrix derivation here is simpler than that 
given in reference 25. We also wish to present slightly more detail in the 
algorithm and to give the corresponding FORTRAN program for an example 
problem. 

As explained in reference 25, a modified treatment of the upper boundary 
results in a regular block coefficient matrix with all the block-diagonal 
elements the same (rather than being different as in eq. (31b)). The modi- 
fied treatment changes the last column of the matrix from the form obtained 
in equation (31b) (equivalent to a Neumann condition applied to Poisson’s 
equation in a certain way, as done in refs. 9 and 17) to another form that is 
equivalent to a Neumann condition applied to Poisson’s equation in a different 
way (as done in ref. 5), which is simpler to treat by cyclic reduction. The 
modification of the upper boundary to accomplish this was originally derived 
by working backward from the desired result. Therefore, the motivation is not 
clear at the beginning of the derivation but becomes clear when a crucial step 
in the algorithm development is reached. 

Consider the mesh shown in figure 2(a), on which the value of m is 
arbitrary and on which the number of y values both for the continuity equa- 
tion and for the rotationality equation is = 2^, where the exponent L 

is an integer. Note that this is one greater than the dimension n used in 
the previous section where figure 2(b) was used for the illustration. The 
rotationality equation is to be applied at the points where u is specified 
on the upper boundary, as well as all the points indicated by crosses. With 
3 defined by equation (25), the assumption of zero values for o& and a“, 
and with the definitions in equations (4) to (8), equations (3a) and (3b) 
become 


For j = 1 to m jand k ~ 1 to n\ : 


g (u . 7 - u. 7 ,)+ v . i - v . 7 . « h s . 7 

0 9 k j-i ,k o ,k y j ,k 

For Q = 1 to m and k = 1 to n: 


U 3 ,k + U J>k+i + '3 ,k V g+\,k ) 


(33a) 


(33b) 


However, for the rotationality equation at the upper boundary only ( k = n^), 
the difference expression for (du/dy)j ^ given in equation (8b) is replaced 
by the first-order-accurate expression: 


(Zu/dy) . - 

c7 9 K 


(U j,k '+ 1 



(34) 
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where, as indicated in figure 2(a), Uj is evaluated on y = y u . Thus, 

as before, ignoring the distinction between j and j 1 or k and for 

indexing, we write 

For j = 1 to 777 and k = ni ; 


'“.I.* + “f.fc+l + 2 


'i-ui 


) = - 


7T h 0) . 7 

2 y dyk 


(33c) 


Although equations (33a) and (33b) are second-order-accurate representations 
of equations (1) at all points, and equation (33c) is only first-order-accurate 
at the upper boundary, the applications for which the algorithm is intended 
are not significantly affected by the lower accuracy at the upper boundary, 
which is assumed to be sufficiently far from regions where any large gradients 
of u and V occur. The factor 1/2 in the denominator in equation (34), 
with the resulting factor of 1/2 multiplying g in equation (33c), is 
crucial for obtaining the desired result. With all specified boundary values 
of Uj ^ an< 3 Vj £, as indicated in figure 2(a), moved to the right sides of 

equations (33a), (33b), and (33c), denote the resulting right side of (33a) 
by fj 9 k (i = 1 to m, k - 1 to n\) and the resulting right sides of equa- 
tions (33b) and (33c) by g • v (j = 1 to m, k = 1 to n x ) . Then with the 
definitions in the section Matrix Definitions and Relations, and with 


3, = -7 = 6 (except B = 0 ) 

1 > d id ^ , 1 

and 


P = q = «i 


(35) 


(36) 


the complete set of finite-difference equations (33) for all j , k can be 
written as 


n\ blocks n\ blocks 

A ■ v / ■ * 


A ,10 





A ! -I I 


u 2 


f 2 

1 

* • * O 


• 


• 

• 1 • • 


• 


• 

A . -I I 


u 


f 

1 


n l 


» i 



— — 


— — 

-II , B 


Vi 


*1 

0 -I ! 

• 1 • 


V 2 


§2 

I B 


\ 


• 

♦ 

r* ““ 1 
•hV 

o 

1 




_ S «i_ 


(37) 
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or, more concisely 


DjU - E T V = f (38a) 

EU + D' V = g (38b) 

D 


Equations (38) can be reduced by eliminating U. Premultiply equa- 
tion (38a) by -E, premultiply equation (38b) by D n ^(A), add the two result- 
ing matrix equations, and use the relations (24a) and (24c) to obtain 



(-1,0, -I) 


+ D f 

n , 


(c)]v = 


(39a) 


where F is defined as in equations (32); but, in this case, p = n\. Again, 
since p = q 9 fp+i must be defined to be zero for equation (32c) to be valid. 


Equation (39a) is a block-tridiagonal equation for V. 
an expanded form: 




m m 


» — 

C -I 


Vi 


Fl 

-I 


V 2 


f 2 

C -I 


• 


. 

-I % c 


V 


F 

m m 


_ «1_ 


. Wl . 


If we write it in 


(39b) 


we can now see the motivation for using equation (33c) containing the factor 
1/2 at the upper boundary. By multiplying the last column of /n- dimensional 
matrices in the block coefficient matrix in equation (39b) by 2 and at the 
same time multiplying the last vector element, V ni , by 1/2, we obtain 


” C -I 


'Vi ' 


“Fi" 

- 1 c 


V 2 


f 2 

-I •. -I 


* 

— 




* 


* 

C -21 


V n 


F n 

-I C 


l*V.„ 


F 

J 1 


“ n \ J 


(40) 


The block-tridiagonal coefficient matrix now has a form for which the cyclic 
reduction is relatively simple for Yi\ - 2^, L being an integer. The cyclic 
reduction of equation (40) to obtain all values of Vj ^ is adequately 
described in reference 25 (pp. 68-72). The values of Uj ^ are then obtained 
directly from the lower half of equation (37) (i.e., the rotationality equa- 
tion for each j ,k) , from which 
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(41a) 


and 


U 

»1 



8 


n l 


U, = U. 


fc +1 


+ B \- 


ft = n to 1) 


(41b) 


where (see eq. (21b)) the j component of the vector represented by each prod- 
uct BVj, is just 3 - t>j+i,fe)» where V m \ 9 k is set equal to zero for 

all k = 1 to rii (since it is already included in 0j 9 k)' 

A FORTRAN program to solve an illustrative example problem, using sub- 
routines corresponding to version A, is described in a later section. 


VERSIONS B AND C — FIRST REVISED AND EXTENDED C AU CHY-RI EMANN SOLVER 


The motivation for developing version B was the prospect of increased 
efficiency over version A. Because version A (as described above) contains a 
block-tridiagonal matrix equation to solve that is equivalent to a discrete 
Poisson equation with a Neumann condition, the cyclic reduction is less effi- 
cient than would be the cyclic reduction of the corresponding equivalent 
matrix equation derivable from Poisson's equation with all Dirichlet condi- 
tions. The latter is known to require one less cycle of reduction containing 
a significant number of operations, in both the forward and backward recur- 
sions. Therefore, the objective was to change the condition application at 
the upper boundary so that the block-tridiagonal matrix obtained would have 
all diagonal blocks the same, as in equation (40), but also have all super- 
and subdiagonal blocks the same, as in equation (31b). The revised method for 
treating the upper boundary could also alleviate effects of the decreased 
accuracy at the upper boundary in version A in problems where those effects 
could be significant. 

At the time the revised version B was being developed to increase the 
efficiency, it was also found, in the studies that were to be reported later 
in references 36 and 37, that certain terms needed to be added to the solver 
to stabilize the iterations in the fast semidirect iterative method for the 
nonlinear, mixed elliptic-hyperbolic problem of transonic flow. Therefore, 
the algorithm for version B was extended to include the terms a + u"*" + 
in equations (3), with constant values of a + and a”. This extended algo- 
rithm, version C, is the one used for all the calculations made in refer- 
ences 36 to 38. Since version B can be regarded as a special case of 
version C, the derivations are combined. However, the subroutines for 
version B are slightly more efficient than those for C and can be used for 
the same problems as version A, where the extra terms are not required. 

Consider figure 2(b), on which m is arbitrary and n\ = n 4- 1 = 2^ 
with L an integer for the simplest and most efficient cyclic reduction of 
resulting equations. Note that on the upper boundary, at k T = n\ (y - y u ) , 
values of u (circled symbols) are specified. It is also indicated that. 
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at a half interval below the upper boundary, at k = n 9 values of V are 
also to be specified. (Such conditions in most physical problems are not 
difficult to obtain. For example, in refs. 36 to 38, analytical asymptotic 
conditions valid sufficiently far from y = y ^ were specified for u at 
k 1 = ni and for V at k = n. In the event that values of u and V would 
be known only at a single boundary line, y - y U9 the discretized continuity 
equation could be used to obtain corresponding values of v on k = n to 
specify as input.) As was the case in the previous section for version A, 
the motivation for this particular modification of the upper-boundary treat- 
ment is not obvious at the beginning since the necessary modification was 
originally obtained by working backward from a desired result. After a cer- 
tain point is reached in the derivation, the reason for specifying V at 
k = n will become clear. 

The specification of V on k = n, the crucial feature of this algorithm 
that leads to the desired result, allows us to add an unknown, Vj > n9 to the 
left side of each continuity equation and to add the corresponding specified 
quantity, denoted by (prp)j 9 to the right side. Thus, even though we specify 
(Vrp)j , we also shall determine v j 9 n as P art °f solution algorithm. Also 

for use in equations (3) with nonzero values of a + and a”, let 

5 1 = ~ h x a+ » “2 E ~ h x a ~ ( 42 ) 

and in terms of these parameters, along with S defined by equation (25), 
let 

hj - ^ - 5 1 > 

B . = 3(1 + a ) (except 3 2 i E 0) 

Then corresponding to figure 2(b), with equation (3a r ) applied at all the dots 
and equation (3b) applied at all the crosses, with the definitions in equa- 
tions (4) to (8), we write for all j = 1 to m and all k = 1 to n: 



ai)u j,k 3(i + a 2 )^._i + 


V 3,k-\ 


+ v . = h s\ 7 + (y_) . (44a) 

3 ,n y 3, k T 3 


-u . « + 


(44b) 


With an appropriate expression for sj j, 9 as in equation (3a), these equations 
are second-order-accurate representations of equations (1) . 


Now, with all specified boundary values of Uj ? ^ and ^j 9 k as indicated 
in figure 2(b) (except Vj n ) moved to the right sides of equations (44), 
denote the resulting right sides, respectively, by fj y and gj ^ (j = 1 to m 9 
k = 1 to n). Then with the definitions in the section ’Matrix Definitions and 
Relations, and with 3i i and So a defined by equations (43) and 

1 9d 9d 

p = q = n (45) 
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the complete set of finite-difference equations (44) for all j 9 k can be 
written as 


n blocks n blocks 



(46) 


With the additional use in this section of the matrix H n defined by equa- 
tion (16) , the more concise formulation of equation (46) is 

D^U + ( -E T + H^)V = f (47a) 

EU + D V = g (47b) 

B 

Equations (47) can be reduced by eliminating U. Premultiply equa- 
tion (47a) by -E 5 premultiply equation (47b) by (A ) 9 add the two result- 
ing equations, and use the relations (24a) and (24b) to obtain 

[T (-1 ,C - 1) - J (I) - EH ]V = F (48) 

n n n 

where F is defined as in equations (32) with p = n 9 and where, since 
q = p, fp+i in equation (32c) must be defined to be zero. At this point, 
the motivation for introducing Vj n into equation (44a) , which resulted in 
the term with the matrix H. n in equation (47a), becomes evident. The iden- 
tity 

t E M (I) ]H n = -J n O) («) 
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yields, from equation (48) 


[T»(~I*C.-I>]V = F (50a) 

which is the desired simple and very regular form of the block-tridiagonal 
equation for V: 



(50b) 


Version C uses cyclic reduction of equation (50b) , with the matrix C 
defined by equations ( 21 c) and ( 21 d) , where &i 9 j anc * are defined by 

equations (43). Version B is for the special case where 04 and a 2 in equa- 
tions (43) are zero. The cyclic reduction process for obtaining all values of 
Vj from equation (50b) is slightly different from that used in version A 

anci is described in appendix F. The values of u j 9 k are then obtained 
directly from the lower half of equation (46), from which 

\ = V k + 1 + " &k (k “ n > n - i, . . 1) (51) 

where ^n +1 is set to zero, and where the j component of the vector 

represented by each product BV^ is just 9 k ~ v j+l 9 k) with 'Om\ 9 k set 

equal to zero for all k = 1 to n (since it is already included in 9j 9 k)' 

FORTRAN programs to solve the illustrative example problem, using sub- 
routines corresponding to both versions B and C, are described in a later 
section. 


VERSIONS D AND E - SECOND REVISED AND EXTENDED CAUCHY-RIEMANN SOLVER 


For current and future applications of a fast direct Cauchy-Riemann 
solver, it is desired to apply conditions on V on the horizontal centerline 
of the mesh (by methods beyond the scope of this report). It is therefore 
desired to have the upper and lower boundaries symmetrically located above 
and below the centerline and to have the same type of boundary condition 
applied at both boundaries. We therefore consider the mesh configuration 
shown in figure 2(c), on which boundary values of V are specified on the 
upper and lower boundaries, u is specified on the left, and v on the right. 

There is a further motivation for pursuing the formulation for the con- 
figuration in figure 2(c). One quickly observes that each column of dots 
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contains one more point for the continuity equation to be solved than the 
number of crosses in a column of the points for the rotationality equation. 
Correspondingly, one observes that the number, q , of unknown values of u 
in a vertical column is also one more than the number, p, of unknown values 
of V . Thus, in equations (11), 


q = p + 1 


(52) 


If one then recalls the formulations in previous sections (e,g., eq. (29)), 
he notes that the quarter partitions of the block coefficient matrix are 
square because q = p there. Consider that the number of rows of mxm 
matrices in the upper partition of equation (29) is q (the number of 
vectors), whereas the number of rows in the lower partition is p (the number 
of Vfc vectors); the number of columns in the left partition of the coeffi- 
cient matrix must be q , whereas the number of columns in the right partition 
must be p in order for the indicated matrix multiplication to be defined. 
Thus the upper left partition has block dimension q x q 9 the upper right qxp , 
the lower left pxp, and the lower right pxp. Recalling also the identity 
(24d) , we have an indication that the simplest and most regular coefficient 
matrix after reduction may be obtained quite naturally for q — p + 1, 
whereas it was not obtained naturally (i.e., without the tricks used for 
versions A, B, C) , for example, in the section Simplest Formulation of Matrix 
Equations, in which it was assumed that the number of rotationality equations 
was the same as the number of continuity equations (p = q ) , and all block 
matrices other than vectors were square. Thus, the simplicity of the matrix 
in the identity (24d) (with the implied corresponding high efficiency of the 
algorithm that would result) encourages the pursuit of the formulation cor- 
responding to figure 2(c) with nonsquare matrices. 

Versions D and E of the extended Cauchy-Riemann solver are derived here 
as a single algorithm for the mesh in figure 2(c). Version D is the special 
case with constant coefficients, a*~ and a”, in the extra terms (a + u + + a~u“) 
in equation (3a), whereas version E allows coefficients that are arbitrarily 
variable in x . 

Figure 2(c) requires 


p=n, q = ni = n + 1 (53) 

This requirement, with n\ = 2^ (L, an integer), will result in a matrix 
equation for which the cyclic reduction is simplest and most efficient. The 
value of m, as in the other versions, is arbitrary. 

For use in equations (3), let a + and a~ be defined, respectively, at 
the same points as i& and vT in figure 1(a). To use the indexing notation 
of figure 1(b), then, let 


= -h a + 
x 


a„ . t = - ha 


(54) 
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and in terms of these parameters, along with 6 defined by equation (25), let 


6 i,j ‘ 6(1 - 


B 2 ,j 6(1 + °2,3-l ) 


(except 



(55) 


Then corresponding to figure 2(c), with equation (3a f ) applied at all dots 
and equation (3b) applied at all crosses, with the definitions in equa- 
tions (4) to ( 8 ), we write 


For all j = 1 to m and all k = 1 to : 


6(1 - a_ .)u . 7 - 6(1 + a. . )u . 7 + V . 7 - v . 7 - h s\ 7 (56a) 

J ,k 2,j-i j-i y k j 9 k j 9 k - 1 y o y k 

For all j = 1 to w and all k = 1 to n: 


-V,* + “j.fcH + 6 <, V,* " Vl.*’ ■ 'W.* 


(56b) 


With an appropriate representation of sj ^ £ as in equation (3a), these equa- 
tions are second-order-accurate representations of equations ( 1 ) . 


With the specified boundary values of u j 9 k and v j,k indicated in 
figure 2(c) moved to the right sides of equations (56), denote the resulting 
right sides, respectively, by fj 9 k (j = 1 to m 9 k = 1 to n\) and gj 9 k 
(j = 1 to m and k =* 1 to n) . Then, with the definitions in the section 
Matrix Definitions and Relations, and with 6 } j and 62 j defined by equa- 
tion (55), the complete set of finite-difference equations (56) for all 
values of j 9 k indicated can be written as 
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or, more concisely, 


D 4 U - E2V = f (58a) 

EU + D g V = g (58b) 

m 

Note that, in this formulation, E and E^ are nonsquare and and are 

square but of different block dimensions (see eqs. (22)). 

Equations (58) can be reduced by eliminating U. Premultiply equa- 
tion (58a) by -E, premultiply equation (58b) by D n (A), add the two resulting 
equations, and use equations (24a) and (24d) to obtain directly the simple 
and very regular form of the block-tridiagonal equation for V: 

[T n (-I,C,-I)]V = F (59) 

where F is defined as in equations (32) with p = n, q - ri\. (In this 
case, there is no special condition on fp+i as there was in the cases where 
E was square.) The expanded form of equation (59) is identical to equa- 
tion (50b) , except that in this case C has been allowed to have variable 
elements (eqs. (21c) and (21d) with (55)). 

Version E of the solver uses cyclic reduction of equation (59). Ver- 
sion D is for the special case where a i 9 n_ anc * a 2 j-l in equation (55) are 
constants. The cyclic reduction process ror obtaining all values of Vj ^ 
from equation (59) is described in appendix F. 

The procedure for obtaining the after all V^, are known is somewhat 

different from the previous versions since k ranges from 1 to n for 
but from 1 to n\ for U^. After U is determined, the lower partition of 

equation (57) (i.e., rotationality equations) can be used to obtain the 

remaining U^: 


U k = U £ +1 + BV /, ~ (k = n, n - (60) 

in which the j component of the vector represented by each product BV^ 
is 3 (Vj ? with v m^ y k set ec I ua ^- to zero for all (since it is 

already included in gj * jf ) . To obtain U n ^ before equation (60) is used, 
several different approaches can be taken. Let us refer to these as options 
(a), (b), (c), and (d) and consider them in detail. For option (a) to 
determine U n ^, we use the bottom row of the upper partition in equation (57), 
that is, the continuity equation for the top row of dots in figure 2(c): 


AU 

n 


1 


- V 

n 


f 

n l 


(61a) 
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or 


u . 

3**1 


U M. + V . + f . ^ 

j V 2,3 3 -l,n l 3 ,n 3 ,n l ) 


(61b) 


for j = 1 to m (with set equal to zero since it has been included in 

„ ) . Because some difficulty has been encountered and studied in use of 

this option in version D of the solver, it has not been included in version E 
with variable 3l,j and 32, j- Therefore, for discussing the nature of the 
difficulty, we consider only the case (i.e., version D only) where 


$1 j 3(1 ot i ) 


3 2 . = 3 2 = 3(1 + a 2 ) (except 0 2 1 
^ y a 9 


0) 


(62) 


with and 6 2 as constants. Further, let 


W . — u. . 

3 3 >^1 


r . = v . + f. 

3 3,n 3 ,n l ) 


and write equation (61b) as 


(63) 


3 1 W j ~ &2 W n-l = r i ^ = 1 to m) 


(64a) 


in which = 0 or, equivalently, as 


w . - w . , = (X - 1 )w . , + — r . U = 1 to m) 
J J-l J " 1 Si J 


where 


(64b) 


A e 3 2 /3 x 


(65) 


This is a common form of difference equation, which can be studied by standard 
techniques in stability theory to determine that the equation is stable if 
| A | < 1, i.e., | S 2 / 3 i | ^ 1. In fact, the difficulty encountered was traced 
to an instability in the determination of the u j 9 rii by this option only 

when the values of 3i and S 2 were chosen so that |s 2 /3i| > 1. 

Because of the stability difficulty with option (a) in these cases, 
further options were considered that find the U- „ without using the top 

row of continuity equations on the mesh in figure 2(c). Option (b) is to just 
specify each Uj , in the step of the algorithm before equation (60) is 

used, as 
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I 


U . = (u m ) ■ 

j,n x T'q 


(j = 1 to w) 


( 66 ) 


where (u%i)j is obtained from some analytical representation such as an asymp- 
totic boundary condition evaluated at each j and at k 1 - Yi\ in figure 2(c), 
Options (c) and (d) use, instead, different variants of the rotationality 
equation to replace equation (56b) for application at points on y = y u 
(fig. 2(c)) halfway between the V input points and directly above each of 
the interior u points (i.e., at k - , all j) . For these two options, 

let (uy be the respective values of u specified at these points on 
y - y U9 and denote du/dy at these same points by (du/dy)^ . Option (c) 
will use a first-order-accurate expression for (du/dy)^ j to replace the 
difference approximation given in equation (8b) that resulted in the rota- 
tionality difference equation (56b). Similarly, option (d) will use a dif- 
ferent second-order-accurate approximation for (du/dy) ^ ^ . The first-order- 
accurate expression, to replace (du/dy)j ^ given by equation (8b), at 
k - ri \ 9 is ’ 

(§),., - w, - /m 

Thus, for option (c) , the replacement for equation (56b) rearranges to 

u . = (u m ) . + 7 T &(v . - v . , _ \ “ h a) . (j = 1 to m) (68) 

T j 2 V j ,n x 3+1 9 nJ 2 y j 9 n x 

For option (d) , the second-order-accurate expression, to replace (du/dy) j ^ 
given by equation (8b), at k = n\ 9 is 


(*),., * [ 8< Vj - + v 

The resulting replacement for equation (56b) is 

8(u~) . - 9 u . + u . - 36 (v . - V . | = - 

T 3 3 ,ni 3 ,n \f 


■3h co . 

y 3 > n i 


(69) 


(70) 


Before equation (70) can be solved for u j t ni> one must eliminate u y ,n’ 
This can be done by adding equation (70) to (56b) written at k = n: 

-u . + u . +3(u. - v . , ) = ~h co . 

3,n 3 >ni 3 ,n 3+1 ,n' y 3 ,n 

to obtain, after rearranging, 

u ‘ ~ + 8 ^(? 3 ‘,n 1 ~ U j+i,n 1 ) + 8 " U j+i ,n^ 




4 (3co . + co . J (3 = 1 to m) 

8 y\ j.ni J,n/ 


(71) 


which is the equation used to evaluate u j 9 n\ f° r °P t i° n (d) before equa- 
tion (60) is used in the algorithm. 
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FORTRAN programs to solve the illustrative example problem, using sub- 
routines corresponding to both versions D and E, are described later. 


EXAMPLE PROBLEM FOR ILLUSTRATION 


To illustrate the use of the variants and extensions of the Cauchy- 
Riemann solver, and as a basis for a complete sample program to be given for 
each version, consider the same problem used in reference 25. The sample 
programs can be used to determine such factors as computing times required 
for significant portions of the calculations and accuracy of the numerical 
results . 

The problem is the small-perturbation formulation for steady, irrota- 
tional, incompressible, inviscid flow over a thin symmetrical parabolic-arc 
biconvex airfoil aligned with a uniform stream far from the airfoil. This 
problem has a simple mathematical formulation and the analytical solution 
is available for comparison with numerical results. Also, the problem has 
the interesting property of containing singularities, which should be captured 
by the numerical solution. The problem formulation and programs also consti- 
tute a basis for extensions to the more complex problems of nonlinear subsonic 
and transonic aerodynamics (refs. 35 to 38). 

In Cartesian coordinates (x 9 y), let U and V be the respective com- 
ponents of velocity. For the classical small-perturbation approximation, let 

V = U (1 + tu) , V = V TV (72) 

oo 7 oo ' 


where u and V are then dimensionless scaled perturbation velocities, U 
is the uniform velocity at infinity, and r is the small thickness ratio of 
the airfoil. If these expressions are substituted into the governing gas- 
dynamic conservation equations, the resulting equations for u(x 9 y) and 
V(x 9 y) are the Cauchy-Riemann equations: 


dU ( dV 

dx 9 y 


0 


(73a) 


9 u 9t? 

9 y 9x 


(73b) 


Let both x and y be normalized by the chord length of the airfoil and let 
the x axis be along the chord and the y axis be the perpendicular bisector 
of the chord. Then the biconvex-airfoil surface y - yjy(x) is given by 

y b (x) = ±t(| - 2x 2 ^ , (- j < x < -|) (74) 


In thin-airfoil theory, the condition of flow tangency at the airfoil surface 
y - z Jfo(x) is transferred to y - 0 by use of Taylor’s series (see, e.g., 
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ref. 43).- The first-order, thin-airfoil result for the boundary condition is 
then 


( 1 1 
- ~2 < x < ~2 

or 

v(x 9 0 + ) = -bx in (“ \ ^ x ^ f) 

= 0 for j rc j > 

The condition of uniform flow at infinity requires that 

u 9 v 0 as x 2 + y 2 -> 00 


(75a) 

(75b) 


(76) 


The analytical solution to equations (73) with conditions (75) and (76) (see 
table A. 2 of ref. 44, p. 21), in terms of the complex variable z = x + iy , 
is 


u - iv 




(77) 


from which, at y - 0, 


u(x 9 0 ) 



x £n 


x + d/2) r 

x - (1/2) L 


(78) 


Note from conditions (75) that V has a discontinuous jump at both the leading 
and trailing edges (|^| = 1/2) and, from equation (78), that u goes to 
negative infinity there. 


For the computational problem, because of symmetry we consider only the 
half-plane y > 0 and use conditions (75) at points along the computational 
boundary y - y^ = 0. For the outer boundaries, we let y = y u , x = Xg 9 and 
x = x u define the boundary lines where conditions are to be applied and con- 
sider using either u = u 0 (x,y) or V = v 0 (x 9 y) as the conditions to be 
applied there, where u Q and V Q are the analytical expressions obtained from 
the real and imaginary parts of the exact solution (77). With these exact 
solutions applied on the outer boundaries, replacing the asymptotic condi- 
tions (76), the illustrations will not be affected by errors due to approxi- 
mate methods for applying conditions at infinity. However, one could also use 
zero perturbations (u = V = 0) on the outer boundaries as approximate condi- 
tions, and this option is allowed for in the programs described in the next 
section. Therefore, let us denote the outer boundaries (x^ 9 x U9 y u ) as B and 
denote the conditions there as 


u = Ug(x 9 y) or v = V^(x 9 y) on B 


(79) 
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where Ug may be either u Q or 0 and Vg may be either V Q or 0. The 

choice of whether to specify Ug or Vg at points on a section of B 

depends on which version of the algorithm described previously is used. 
Furthermore, in versions B and C, where (.Vrg)j must be specified at a half- 
interval inside the upper boundary, (Uy)j* can be obtained from the function 
Vg(x,y); in versions D and E, in which (wy)j' must be specified either at a 
half-interval inside the upper boundary for option (b) or on the upper 
boundary for options (c) and (d) , (uy)^* can be obtained from the function 
u B (x,y ). 

To test and illustrate all versions of the algorithm, including those 
with the extra terms in equation (3a) , using this example problem governed by 
equations (73), let us represent equations (73) by the finite-difference equa- 
tions (3) with the notation defined in equations (4) to (8) and with s and w 
equal to zero, so that, in equations (3a) and (3a'): 

s’ (•) = a + Up + cTu 2 (80) 

in those cases for which a + and a~ are not zero. Therefore, for version C, 
the term Sj ^ in equation (44a) is 

h s'. , = -&a 2 (u ). - 8a \(u ) . (81a) 

y 3,k ° J-1 ,k o a ,k 

and for versions D and E, in equation (56a), 


h s'. 7 = -8a . (u ) . - 8a . ( u ). 

y 0 ,k 2,3-1 o j- i,fc 1,3 


where u Q is the exact solution defined above. 


(81b) 


To summarize the numerical example problem for illustration, the differ- 
ence equations approximating the partial differential equations (73) are 
equations (3) with (80) represented by equations (81) for the algorithm ver- 
sions with the extra terms; the boundary conditions are the discretized 
applications of equations (75) on y = y % = 0 and of equations (79) on 
x ~ x = x U y y = y u , with the choice of whether u or V is specified 
on any section of B depending on the algorithm version being used. 


After the example problem is solved by any version of the algorithm, it 
is of interest to determine u on y = 0 (which approximates the airfoil- 
surface velocity for < x < ^ in thin-airfoil theory). Since u on 
y - 0 is not determined directly in the solutions obtained on the meshes in 
figure 2, a second-order-accurate representation of the rotationality equation 
can be applied on k = 0, with one-sided y differencing, analogous to equa- 
tion (70) used in versions D and E for k = The rotationality equation 

can be arranged to obtain (ref. 25): 


"J,0 ■ (1/8)[9t V, 1 


- U 


3,2 


36(t V,o 


- V . ) 

3+ 1,0' 


3 h a) . ] 

y 3,o J 


( 3 = 1 to m) (82) 
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(in which u j 9 o is the value at k = 0 rather than at fc 1 = 0). In this 
example problem, of course, uy = 0 in equation (82). 

FORTRAN PROGRAMS AND RESULTS FOR THE EXAMPLE PROBLEM 
USING VERSIONS A TO E OF CAUCHY-RIEMANN SOLVER 


The purposes of this section are (a) to describe FORTRAN programs, listed 
in appendices A to E, for solving the simple example problem outlined in the 
above section, using versions A to E of the Cauchy-Riemann solver and (b) to 
describe results of calculations for the example problem. 


FORTRAN Program Descriptions 

The FORTRAN programs are written for use on a Control Data 7600 com- 
puter. Slight modifications would be required for use on most other com- 
puters that have FORTRAN compilers. Each program consists of a main program 
and three subroutines, as listed in table II. The last letter of the name of 


TABLE II.- PROGRAMS AND SUBROUTINES 


UBIC1A 

UBIC1B 

UBIC1C 

UBIC1D 

UBIC1E 

CAL CD 1 

CALCD 

CALCD 

CALCD 

CALCD 

CR2N1 

CR2DUV 

CR2UVC 

CR2UVC 

CR2UVE 

GE2N1 

GE2N1 

GE2N2 

GE2N2 

GEUVE 


the main program denotes the version of the Cauchy-Riemann solver used; thus, 
UBIC1A is the main program to solve the example problem for the biconvex air- 
foil using version A of the solver. The subroutine CALCD used with versions B 
to E computes the quantities d £ m ( ec l s * (F7) in appendix F) . The correspond- 
ing computation in version A is done by subroutine CALCD1. The subroutines 
CR2DUV , CR2UVC, and CR2UVE perform the cyclic reduction as outlined in appen- 
dix F for versions B, C, D, and E. The subroutine CR2N1 performs the version 
of cyclic reduction described in reference 25 for version A. Within each of 
the subroutines for cyclic reduction, another subroutine is called to solve 
tridiagonal equations using the Thomas algorithm (e.g., see ref. 45) for 
Gaussian elimination. The slightly different subroutine versions for this 
algorithm are named GE2N1 , GE2N2 , and GEUVE. 

The use of these partially modular forms of the programs, with separate 
subroutines to perform major steps in the program, is not the most efficient 
but is best for understanding, modifying, and adapting the programs. 

Timer - Each main program includes the use of a timing subprogram from the 
CDC-7600 program library that returns the CPU time T in seconds (from the 
start of the job) either by the statement CALL SECOND (T) or, for example, by 
CPU=SECOND (T) . This subprogram is used to measure three times in each main 
program: T1 is the time for all preliminary computations, excluding input. 
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before the cyclic reduction; T2 is the time for the cyclic reduction to be 
completed (determination of V after F is known) ; and T3 is the time 
required to obtain U after V is known. When the program is used on other 
computers, all statements containing the FORTRAN variables TT1, TT2, Tl, T2, 
and T3 should be removed or appropriately replaced. 

FORTRAN parameters and variables - Before further description of the pro- 
grams, it is convenient to list in table III some FORTRAN parameters and 
variables that either are used frequently or are special and that correspond 
to the indicated equivalent algebraic quantities. For this we note that, 
since FORTRAN does not permit the index zero in an array, the FORTRAN indices 
J and K, for example, are shifted by one from the values of j and k. The 
variable J may represent either JU or JV in table III and the variable K 
may represent either KU or KV in the same sense as j and k are used above 
as indices to represent, respectively, either j or j f and either k' or k. 

Statement functions - At the beginning of each main program, statement 
functions are included to determine X = XJU(J) or XJV(J) and Y = YKU(K) or 
YKV(K) from equations (4) (see table III). Statement functions are also 
given for the exact analytical solutions for u and v 9 UEX(X,Y) and VEX(X,Y) , 
which include within the statements the additional statement functions 
RHO(X,Y) and ANPHI(X,Y). The functions UEX and VEX represent the real and 
imaginary parts of equation (77). In addition, version E contains the state- 
ment functions Fl(JU) and F2(JU), which are arbitrarily specified functions 
to represent a l and a 2 9 j> respectively, in the test calculation. 

Computation of boundaries - Because of the staggered mesh, in some cases 
it is most convenient to specify nominal values of x^ , x u , and y U9 and then 
compute final values of those parameters for which convenient mesh intervals, 
h x and hy , are obtained. For this reason, after the nominal values are input, 
and mi and n\ are input, h x and hy are obtained from equations (9d) and 
(9e) , then equations (9a) and (9b) or (9c) are used to determine the final 
x u> anc ^ Mu* taking into account the fact that it is desired to have 
Xj - -0.5 (as given by eq. (4a)) coincide with the leading edge of the airfoil 
at some J = JLE. This computation is done just after the input is read at 
the beginning of each main program. 

Input - As indicated near the beginning of each main program, versions A, 

B, and E each require only one input card per case (READ statement labeled 2 
with corresponding FORMAT statement labeled 1), and versions C and D require 
two input cards per case (READ statements labeled 2 and 104 with corresponding 
FORMAT statements labeled 1 and 103). In each program, the first (or one only) 
input card format for each case requires the following specifications: 

NEX Integer in column 5. Specification is arbitrary and superfluous unless 
NX is specified as some integral power of 2 (NX=2**NEX) . In the latter 
event, specify NEX equal to log 2 NX. 

NX Integer, right-justified at column 10: Mesh number for x direction, 

equal to 777 * . Specify integer ranging from 4 to 128. (Recommended 
values either have a factor of 10 or are a power of 2; see table I.) 



TABLE III.- FORTRAN SYMBOL DEFINITIONS 


FORTRAN 

Algebraic 

Reference equation 

Reference page 

JU 

0 + i 

(4a) 

9 

JV 

+ 1 

(4b) 

9 

KU 

k.' + 1 

(4b) 

9 

KV 

k + 1 

(4a) 

9 

NX 

mi, m + 1 


8 

NY 

rii, n + 1 


8,9 

JM 

mi + 1 


8 

KM 

rii + 1 


9 

NEY 

L 


9 

XL, XU 

v x w 


6,9 

YL, YU 

y u 


6,9 

XJU(J) 

x o 

(4a) 

9 

XJV(J) 

x r 

(4b) 

9 

YKU (K) 

y'r 

(4b) 

9 

YKV(K) 

yk 

(4a) 

9 

HX, HY 

V \ 

(9) 


U(J,K) 

u 3' ,k 

(7) 

9 

V (JjK) 

V 3,k' 

(7) 

9 

UEX(X, Y) 

u Q (x y y) 
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VEX (X, Y) 

v Q (.x,y) 


31 

BETA 

6 

(25) 


BETA1 

*l, d 

(43) ,(55) 


BETA2 


(43), (55) 


ALPHA1 

a., , a. . 

(42), (54) 


ALP HA 2 

a , a 
2* 2 , J 

(42), (54) 


A 


(21d) 


A1 

Bi 

(21d) 


ALFA 

^3 

(FlOc) 


ALFA1 

e’i 

(FlOc) 


GAMMA 

-a . 

J 

(21d) 


DELTA 

~ y 3 

(21d) 


ALF(J) 

a . 

j 

(21d) 


BET (J) 


(21d) 


GAM(J) 

y 3 

(21d) 


VT (J) 

(v T ) 


22 

UT (J) 
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NEY Integer in column 15. Since NY is always an integral power of 2 
(NY = 2**NEY), specify NEY equal to log 2 NY. 

NY Integer, right-justified at column 20: Mesh number for y direction, 

equal to n\. Specify integral power of 2 (NY=2**NEY) with NEY ranging 
from 2 to 7. 

On this same card, as described in the previous subsection (see table III for 
correspondence between XL and x% 9 etc.), nominal values are specified for 
XL, XU, YL, and YU. The values have E fields of width 10 and are right- 
justified, respectively, in columns 30, 40, 50, and 60. Recommended values 
for the example problem are, respectively, -1.E0, 1.E0, 0.E0, and 2.E0. (For 
the example problem as formulated, YL is always 0.E0.) The same input card 
in all the program versions also requires specification of the FORTRAN 
parameter NBC in column 65 as follows: 

NBC = 1 results in the outer boundary condition being specified as the exact 
analytical solution; 

= 2 results in zero values being specified for the conditions on the 
outer boundaries. 

In addition to the above parameters, which are all specified on the first 
(or only) input card for each case, the following parameter options apply as 
indicated : 

LASCAS = 0 causes the program to return, after a case has been run, to the 
first READ statement for a new case; 

= 1 causes the program to end; should be specified for last case in a 
job. 

NUBDRY =0 in program D, specifies option (a) (eq. (61b)) for determination 
u j 9 ni* 


= 1 

in 

programs 

D 

and 

E, 

specifies 

option 

(b) 

(eq. 

(66)) 

for 

u d,n 1 ; 

= 2 

in 

programs 

D 

and 

E, 

specifies 

option 

(c) 

(eq. 

(68)) 

for 

u j ,ni ; 

= 3 

in 

programs 

D 

and 

E, 

specifies 

option 

(d) 

(eq. 

(71)) 

for 

u 3 ,n x ' 


ALPHA1 and ALPHA2 are arbitrarily specified numbers in programs C and D except 
that, to keep the difference equations elliptic, require both (-°° < < 1) 

and (-1 < a 2 < 00 ) and, in addition, to keep the calculation stable for the 
option NUBDRY = 0 in program D, require also | (1 + a 2 )/(l - ot j ) | < 1. 

In programs A and B, LASCAS is an integer in column 70; in program E, it 
is in column 75; in programs C and D, it is in column 25 of the second input 
card for each case. The integer parameter NUBDRY is in column 70 of the first 
input card for each case in programs D and E. The parameters ALPHA1 and 
ALPHA2 for programs C and D are specified with E fields of width 10 and are 
right- j us tified , respectively, at columns 10 and 20 of the second input card 
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for each case. For convenience, table IV summarizes the input card and 
column locations for the input parameters. 

TABLE IV.- INPUT CARD AND COLUMN LOCATIONS FOR INPUT PARAMETERS 
(Column of right-justification) 


Input 

parameter: 

NEX 


NEY 

£ 

>< 

S 

►j 

>-< 

& 

>-* 

NBC 

LAS CAS 

NUBDRY 

ALPHA1 

ALPHA2 

Version A 














column 

5 

10 

15 

20 

30 

40 

50 

60 

65 

70 




Version B 














column 

5 

10 

15 

20 

30 

40 

50 

60 

65 

70 




Version C 














card, 

column 

1,5 

1,10 

1,15 

1,20 

1,30 

1,40 

1,50 

1,60 

1,65 

2,25 


2,10 

2,20 

Version D 












* 


card, 
co lumn 

1,5 

1,10 

1,15 

1,20 

1,30 

1,40 

1,50 

1,60 

1,65 

2,25 

1,70 

2,10 

2,20 

Version E 














column 

5 

10 

15 

20 

30 

40 

50 : 

60 

65 

75 

70 




Note that the arbitrarily specified statement functions F1(J) and F2(J) 
determine and a 2,j in program E, so these quantities are not input for 

version E; these specifications are easily changed by simply replacing F1(J) 
and F2(J) in the main program by suitable statement functions. These speci- 
fications should satisfy (-» < 04 • < 1 ) and (-1 < a 2 A < 00 ) at all j in 

(1 < 3 < m). 

Outl'ine of 'programs - The FORTRAN programs have the following major steps 
(see listings in appendices): 

1. Specify DIMENSION and COMMON statements. 

2. Specify statement-function declarations. 

3. Read input as described in above section. 

4. Calculate x £ and x u (and y u in versions B and C) from the nominal 
input values. 
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5. 


Calculate needed parameters. 

6. Calculate in versions B and C (see p. 22). 

7. Initialize interior values of fj 9 k a ^d 9j,k ^ to zero f° r version A; 

to plus the right side of equation (81a) for versions B and C; to the 

right side of equation (81b) for versions D and E) . 

8. Calculate needed boundary values of u and V . 

9. Calculate (urp)j in versions D and E. 

10. Complete the determination of fj 9 k anc * by "modifying the 

fringes” to include the boundary values of u and v (see definitions of fj 9 k 
and gj ^ for each version; e.g. , p. 19). 

11. Set appropriate boundary values of Uj and Vj ^ to zero as 
required by the reduction algorithms. 

12. Determine each Fj ^ (denoted F0 in the program listings) accord- 
ing to equations (32) as described for each program version. 

13. ‘Determine needed values of dn m (eqs. (F7)) by calling either CALCD 
or CAL CD 1 . 


14. Determine and Bj (A and Al; see table III above) in all versions 
except E (in which these quantities are in the array BET(J), included under 
step 5 above), 

15. Perform cyclic reduction according to the algorithm outlined in 
appendix F for versions B to E or in reference 25 for version A (call sub- 
routine names beginning with CR in table II) to obtain all unknown values of 

v 3,k' 

16. Determine all Uj ^ by procedure described previously for each 
version. 

17. Reload boundary values of Uj ^ and Vj ^ that were previously set 
to zero. 

18. Print as output, for appropriate values of j and k : J, K, X and Y 

(for U) , U, X and Y (for V), V, UEX and VEX evaluated at same points as U and 
V, and ERRU, ERRV, which are the differences between U and UEX and between V 
and VEX. 

19. Print Tl, T2, T3. 

20. Print values of u j 9 o calculated from equation (82) along with 
values of v j 9 o anc * the exact solution for comparison. 
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21. For convenience in plotting the exact solution for u on y = 0 
near the leading- and trailing-edge singularities, also calculate and print 
values at small intervals near these edges. 

Results 

The results of computations for the illustrative example problem of the 
biconvex airfoil include, in particular, plots of the computed values of 
U on y = 0 versus x 9 compared with the analytical solution, and the comput- 
ing times . 

For this example problem, there is very little difference in the results 
for u on y = 0 computed by the different program versions, A to E, with 
the same input conditions, and for any of the options (a) to (d) in version E. 
The very small differences are virtually indistinguishable on a plot such as 
shown in figure 3, for which the results were computed with NX = 40, NY = 32, 
NBC = 1, the nominal values of x^ , x u , y ^ 9 and y u as recommended above, and 
with various values of and a 2 in versions C and D. For these input 

parameters, the mesh intervals h x and hy are, respectively, 0.05 and 0.0625. 
The outer computation boundaries are at or about 1/2 chord length upstream 
and downstream from the airfoil edges and 2 chord lengths above the airfoil. 
Because of the discontinuous conditions on V at the airfoil edges, 

( x 9 y ) = (±0.5, 0.0), the analytical perturbation solution for u goes to 
negative infinity there. These singularities are captured well by the numer- 
ical algorithms and, even on the relatively coarse mesh used for figure 3, the 
values very near those points are accurate. 



X 


Figure 3.- Perturbation velocities for thin parabolic-arc biconvex airfoil 
(NX = 40, NY = 32, exact conditions on boundaries). 

The computing times Tl, T2, and T3, as described under the subsection 
"Timer,” are listed in table V for various meshes and for all program ver- 
sions A to E. Note that, when used in an iterative calculation such as 
described in the Introduction and the section that follows it, no significant 
part of the time Tl would be done more than once (at the beginning of 
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the program). Thus, since T3 « T2, the only time of real significance for 
practical applications is T2. As one example, we see that the time T2 for the 
cyclic reduction using version C for NX = 40, NY = 32 is 16 msec. In refer- 
ence 37, it was noted that the direct solver for this case required only 
14 msec. The difference is due to the fact that a special tridiagonal- 
solution subroutine that takes advantage of the vector capabilities of the 
CDC-7600 was used for the quoted slightly lower time. (The time saving was 
much more significant when a previous version of the FTN compiler was being 
used, before the FTN 4.4 version indicated in table V became available.) 


TABLE V.- COMPUTING TIMES FOR BICONVEX AIRFOIL USING VARIANTS OF 
CAUCHY-RIEMANN SOLVER ON CDC-7600 
(FTN 4.4, OPT = 2 Compiler) 


Version of 
program 

NX 

NY 

Tl, sec 

T2, sec 

T3, sec 

A 

10 

8 

0.005 

0.001 

0 


20 

16 

.006 

.005 

0 


40 

32 

.009 

.022 

.001 


60 

64 

.013 

.076 

.004 


100 

128 

.025 

.287 

.011 

B 

10 

8 

.001 

.001 

0 


20 

16 

.002 

.004 

0 


40 

32 

.005 

.015 

.001 


60 

64 

.010 

.056 

.003 


100 

128 

.025 

.218 

.010 

C 

10 

8 

.002 

.001 

0 


20 

16 

.009 

.003 

0 


40 

32 

.032 

.016 

.001 


60 

64 

.092 

.058 

.003 


100 

128 

.299 

.223 

.010 

D 

10 

8 

.002 

.001 

0 


20 

16 

.009 

.003 

.001 


40 

32 

.032 

.016 

.001 


60 

64 

.092 

.057 

.003 


100 

128 

.296 

.223 

.010 

E 

10 

8 

.003 

.001 

0 


20 

16 

.009 

.003 

0 


40 

32 

.032 

.018 

.001 


60 

64 

.092 

.062 

.003 


100 

128 

.299 

.240 

.011 
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User’s Choice of Algorithm Versions 

To aid the prospective user in selecting which version of the algorithm 

to use, we note the following features: 

Version A : 

(1) solves the nonhomogeneous Cauchy-Riemann equations; 

(2) uses a grid that is nonsymmetrical in both x and y (conditions on u 
imposed on left and top, V on bottom and right); 

(3) has a potential source of errors at the upper boundary due to first-order 
accuracy of the difference equation there; 

(4) is the least efficient of the algorithm versions; 

(5) was used for iterative computations in reference 35. 

Version B : 

(1) solves the nonhomogeneous Cauchy-Riemann equations; 

(2) uses a grid that is nonsymmetrical in both x and y (conditions on u 
imposed on left and top, v on bottom and right); 

(3) is the most efficient version for the simplest nonhomogeneous Cauchy- 
Riemann equations with the nonsymmetrical grid. 

Version C : 

(1) solves extended nonhomogeneous Cauchy-Riemann equations, with the extra 
terms having constant coefficients; 

(2) uses a grid that is nonsymmetrical in both x and y (conditions on u 
imposed on left and top, V on bottom and right); 

(3) is the most efficient version with extra terms (constant coefficients) 
and the nonsyminetrical grid; 

(4) was used for iterative computations in references 36 to 38. 

Version D : 

(1) solves extended nonhomogeneous Cauchy-Riemann equations, with the extra 
terms having constant coefficients; 

(2) uses a grid that is symmetrical in y (conditions on u imposed on left, 
V on top, bottom, and right; 

(3) is the most efficient version with the symmetrical grid. 

Version E : 

(1) solves extended nonhomogeneous Cauchy-Riemann equations, with variable 
coefficients in the extra terms; 

(2) uses a grid that is symmetrical in y (conditions on u imposed on left, 
v on top, bottom, and right); 

(3) is the only version presented having variable coefficients in the extra 
terms. 
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CONCLUDING REMARKS 


Revised and extended versions of the original Cauchy- Riemann solver of 
reference 25 have been presented, including detailed derivations of the 
algorithms, as well as descriptions and listings of FORTRAN computer programs 
These programs are believed to be simple and adaptable enough to be useful to 
nonspecialist users. 

The results of the computations for the simple example problem are 
believed to demonstrate significant efficiency, accuracy, and potential 
utility of the algorithms. 

In addition to the previous applications of the original and the first 
revised and extended Cauchy-Riemann solvers in semidirect iterative methods 
for the simplest transonic-flow calculations that have been described in 
references 35 to 38, it is anticipated that the newest version, version E, 
will be useful for more general applications in transonic-flow calculations, 
including lifting airfoils. The Cauchy-Riemann solvers are sufficiently 
general that further significant applications can also be expected in other 
fluid-dynamical problems, as well as in broader areas of mathematical physics 
It is expected that further extensions of these methods in the future, to 
include variable meshes and three dimensions, will also be worthwhile. 

Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif. 94035, December 16, 1976. 
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APPENDIX A 


FORTRAN LISTING OF EXAMPLE PROGRAM AND SUBROUTINES FOR 
ORIGINAL CAUCHY-RIEMANN SOLVER, VERSION A 



C MAIN PROGRAM, IIBIC1A76 FOR CDC 7600 

C+++ SOLUTION OF CAUCHY-R I EMANN EQS . BY CYCLIC REDUCTION FOR INCflM. FLOW 
C OVER THIN BICONVEX AIRFOIL. AIRFOIL CHORD I? FROM X=-.5 TO *.5. 

C 

C THE METHOD AND THE PROGRAM ARE DFSCRI BED IN NASA TN D-7934 

C BY E. D. MARTIN AND H. LOMAX. 

C 

DIMENSION U(129»130) ,V(129,130) 

COMMON NEX,NX,NEY,NY,BETS0,A,A1,ALFA, ALEA1,K,V 
C 

XJUIJU) = XL + HX*FLOAT( JU-1) 

YKUIKU ) = YL + HY*(FL0A T (KU)-1.5) 

XJV(JV) = XL + HX *( FLOAT (JV)-1.5) 

YKVIKVK = YL ♦ HY+FLOAT (KV-l) 

RHO(X,Y) = SQRTl ((X+.5)**2 + Y*Y) /( ( X-.5)**2 + Y*Y ) ) 

AN PH I ( X ,Y) * AT AN2(Y »X-»5) - AT AM 2 ( Y, X+. 5 ) 

UE X( X »Y) = FCTR* ( 1. - X*ALOG(RHO(X,Y) ) - Y*ANPHI(X,Y)> 

VEX (X » Y ) = FCTR*{ Y*AL0G (RHO ( X » Y ) )- X*ANPHI ( X, Y) ) 

C 

FCTR = 4. /3. 14159265 

1 FORMAT (415, 4E10. 0,215) 

2 READ (5,1) NEX,NX,NEY,NY,XL,XU,YL,YU,NBC,LASCAS 
CALL SECOND(TTl) 

C THE INPUT VALUES OF XL AND XU ARE NOMINAL VALUES, USED TO 

C COMPUTE THE P XACT VALUES. 

HY = (YU - YL ) /FLOAT( NY ) 

HX = (XU - XLI/FLOAT ( NX ) 

JLF = 1.005 + (-.5 - XU/HX 
XL = -.5 - HX* ( JL E - 1) 

XU * XL + HX* (NX - .5) 

JM = l + NX 
KM = 1 .+ NY 
KMP 2 = KM + 2 
BETA = HY/HX 
BETSO = RFT A**2 

C++* INI T IALI 2F INTERIOR VALUES OF MASS SOURCES AND VORT TC T T Y (F IS 
C EQUIV. TO V AND G IS EQUIV. TO U)~ 

DO 3 K * 1,1^0 
DO 4 J = 1,129 
VIJ,K) = 0. 

U ( J , K) = 0. 

4 CONTINUE 

3 CONTINUE 

C++++ LOAD BOUNDARY VALUES OF U AND V — 

C IF NBC IS 1, LOAD EXACT (.1 AND V; IF NBC IS 2, LOAD U=0 AND V=0 FAR 

C FROM AIRFOIL. 

K = KM + 1 
DO 5 J = 2, NX 

V(J,1> = 0. 

X = XJV(J) 

IF (X*X.LE.. 25001) V(J,1)= -4.*X 
IF (ABS(X*X-.25 ).LT.l.E-5) V(J,1) = .5*V(J,1) 

IF (MBC.FQ.ll U ( J , K I = IIEX(XJU( JI,YU) 

IF (NBC.EQ.2) U(J,K) = 0. 

5 CONTINUE 

DO 6 K = 2, KM 

IF (NBC.EQ.2) GO TO 33 
U( 1 , K) = UEX ( XL »YKUt K) ) 

V(JM,K) = VEX(XU,YKV(K) ) 
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GO TO 6 

33 IJ(lfK) = 0. 

VCJM,K) = 0. 

6 CONTINUE 

C+++ MODIFY FRINGE? OF INTERIOR to INCLUDE BOUNDARY VALUES — 

DO 7 J = 2, NX 

VI J * 21 = VI J *2 ) + V(J,I> 

U( J » KM I = U(J,KM| - UCJ.KM + II 

7 CONTINUF 

DO 8 K = 2, NY 

V C 2* K I = VC 2» K I + RETA*UC1,K) 

U ( NX »K I = U C NX t K I + BETA*VCJM,K) 

8 CONTINUE 
K = KM 

VC 2 »K ) = V C2 »K I + BET A*U C 1»K I 
UCNX.K) = UC NX »K ) 4- . 5*BETA*V C JM, K l 
C+++ ZERO APPROPRIATE BOUNDARY VALUES — 

K * KM + 1 

DO P J = 1,JM 

UC J * K ) = 0. 

VC J » K> = 0. 

VCJ,1I = 0. 

P CONTINUE 

DO 10 K = 1 * KM 
VCJM,K) = 0. 

UC1 ,K) = 0. 

VC 1 » K) = 0. 

10 CON T IMtJE 

C+++ NEXT DETERMINE VFCTOP FO . 

DO 11 K = 2, KM 
DO 12 J = 2, NX 

VC J * K) = V C J ♦ K ) - V C J » K+l ) + RFTA*CIJCJ,K) - UCJ-1,KI) 

12 CONTINUE 

11 CONTINUE 

CALL CALCD1CNFYI 
A = 2.*C 1. + BFTSOI 
A 1 = 2 . + BETS Q 
' TT 2 = tti 

T 1 = SECGNDCTTI) - TT? 

CALL CR2N1 
T r ?. = tti 

T 2 = SECONDCTtd - T t? 

C+ + + MOW V IS DETERMINED FOR J = 2 TO NX AND K = 2 TO MY + 1. 

C NEX T DETERMINE U p 0P THF SAME J AND K. 

BE TK = .5*BET A 
DO 15 KP = 2, KM 
K = KMP2 - KB 
IF (K.LT.KM) BETK = RETA 
DO 16 J = 2, NX 

U C J»K I = UC J t K+ 1 1 - UCJ.K) + BFTK*CVCJ,K) - VCJ+1,K>) 

16 CONTINUF 

15 CONTINUE 
TT 2 = T T 1 

T3 = SECONDCTTl) - TT 2 
C+++ RELOAD BOUNDAPY VALUER n F U AND V. 

C IB NBC IS 1, LOAD EXACT U AMD V: IF NBC IS 2t LOAD U=0 AND V=0 FAR 

C FROM AIRFOIL. 

K = KM + 1 
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DO 17 J - 2 t NX 
V(J,1) = 0. 

X « XJV( J ) 

IF IX*X.LE*. 25001) V(J,1 ) = -4.*X 

IF <ABS< X*X-.?5).LT.l.E-5) VCJ,1) = .5*V(J,1) 

IF (NBC# EQ.l ) U ( J » K ) = UEX(XJU< J) ,YU) 

TF INBC.EQ.2) U(J,K) = 0. 

IT CONTINUE 

DO 18 K * 2, KM 

IF (NRC.E0.2) GO TO 34 
UfltKI = UEX ( XL »YKU( K ) ) 

V(JM f K) = VEXCXU ? YKV(K) ) 

GP TO IB 

34 UCltKI = 0. 

V(JM,K) = 0. 

IB CONTINUE 

19 FORMAT <1H1*$0LUT ION OF CA UC HY-R I FMANN FQS. BY CYCLIC REDUCTION.*/ 
1* INCOMPRESSIBLE FLOW OVER THIN BICONVEX AIRFOIL. CHORD IS FROM X= 
2-. 5 TO +*5 • *//* NBC=*t I 3>*. (IF NBC IS l, EXACT BOS ARE IMPOSED 
3FAR FROM AIRFOIL; TF NBC IS 2, U AND V APE ZERO ON OUTER BOUNDARY. 
4*// 

5* NX=*,I3,*t MY=*fI3f*t XL = *tF10.6f*t XU=* » FI 0. 6 ,* , YL=*,F10.6, 

6*, YU=* t FI 0 .6 t ** HX=* ?E 15 .7 » * t HY=*, El 5. 7, ***//9X f 1HJ , 3X, 1HK f 6X, 
73H XJ U, 7X,3HYKU,8XtlHU,BX,3HUFX,8X,4HERRU,llX t 3HXJVf 7X t 3HYKVt 8X,1HV 
8 » 8 X ? 3 HVEXt 8X,4HEPRV) 

WRITE! 6, 19) NBC , NX, NY, XL , XU f YL , YU , HX , HY 

20 FORMAT f6X, 2 I 4, FI 1.6, 3*10. 6, FI 5. 7, FI 1 . 6 ,3F1 0. 6 , FI 5 . 7 ) 

KI = NY/ 16 

IF (KI.EO.O) KI ^ 1 
IF (NX-?**NEX.EQ.O) GO TO 21 
JI = NX/20 
GO TO 22 

21 JI = NX/16 

22 IF CJI.EQ.O) JI = 1 

23 rpRMATf IX) 

DO 24 K = 1 , KM, KI 
WRTTE(6,23) 

YA = YKU(K) 

YB = YKVCK) 

DO 25 J = 1 1 JM ♦ J! 

XA = XJU(J ) 

XB = XJV(J) 

IF (J.EQ.JM.OR.K.EQ.l) GO TO 26 
UEXA = UEX { XA , YA ) 

UA * U(J,K) 

GO TO ?7 

26 UEXA = 0. 

UA = 0. 

27 IF (J.EQ.l) GO TO 28 

IF (K.NE.l) GO TO 35 

IF ( A8S ( XB*X B— .25 ) .GE.l.E-5) GO TO 35 
VEX* = 0* 

GO TO 36 

35 V C X5 = VFXCXR,YB) 

36 VB = Vf J * K ) 

GO TO 29 

28 VE XB =0. 

VB = 0. 

29 FRRIJ - tJ A - UEXA 
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ERRV = VB - VEXB 

WRITEl 6,20) J »K , XA , YA ,UA ,UE XA , EPRU, XB, YB » VB , VEX B, FRRV 
25 CONTINUE 

24 CONTINUE 

30 FORMAT ( //* T 1=* , E 12*4, * SEC., T2=*,E12.4,* SEC., T3=*,F12.4,* SEC. 

I* ) 

WRITE ( 6, 30 ) T 1, T 2 ,T3 

38 FORMAT C ///* VELOCITIES AT Y = 0.*) 

39 FORMAT C//26X,IHJ, 6X, 3HX JU , 1 OX, 1HU ,1 3 X , 3HU C X ,11X,3HXJV ? 10X , IHV , I3X, 
13HVEX//) 

40 FORMAT C23X, I4,F11 .6,2E15.7,F11.6,2F15.7) 

WRITE (6,38) 

WRITE (6,39 ) 

DO 42 J = 2, NX 
XA = XJU(J) 

XB = XJV(J) 

UE XA = 0. 

VEXB = 0. 

IF (ARS(XA*XA-.25) .GE.l.F-5) UFXA = UEX(XA,0.) 

IF ( XB*XR.LE.. 25001) VEXB = -4.*XR 
VB = VCJ, 1) 

UA = 0.125*C9 .*U( J,2)-U( J, 3)+3.*RFTA*( VC J, 1 )-V( J+l , 1) ) ) 

WR ITF(6,40) J , XA , !JA , UE XA * XB , VB , VEXB 
42 CONTINUE 

WR T TE ( 6,45 ) 

DELTX * .1*HX 
X = — • 5 - 1.1*HX 
N = 0 

46 N = N + 1 

IP (N.FQ.2) X = .5 - 1 . 1 * HX 
WR I TE ( 6 , 23 ) 

DO 47 I = 1,21 
X = X + DELTX 
UF XI = 0. 

VEX1 = 0. 

IF ( AB SC X*X-.25).GE.l.E— 5) UFXI = UFX(X,0.) 

IF (X*X.LE. .2*^001) VEX1 = -4.*X 
WRITE (6 , 48 ) X , U EX 1 , V FX 1 

47 CONTINUE 

IF (N.EQ.l) GO TO 46 

45 FORMAT ( //* EXACT VELOCITIES AT EDGES.*// 29 X, IHX, 11X,3HU P X, I2X, 3HVE 

IX) 

48 FORMAT (23X, FI0.6, 2F15.7) 

IF CL ASC AS • EO .0 ) GO TO 2 
STOP 

END 


SUBROUTINE CALCDl(NEY) 

C DIMENSIONS OF D (L » M ) ARE LMAX = NEY AND MM A X = 2**LMAX. 

DIMENSION D C 7, 1 28 ) 

COMMON /DC/D 
LMAX * NEY 
DC 1,1) = SORT C 2. ) 

DC 1,2) = -DCl, II 
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10 


20 

30 


o++ 

c 

c 


c 


O'M- 


2 


3 

1 


on 30 L a 2 , LMAX 
MMAX = 2**L 
MMAXH a MM AX/2 
MMAXH1 = MMAXH + 1 
00 10 M = 1, MMAXH 

D ( L , M ) = SQRT<2.+ 0(1-1, M 1 > 
CONTINUE 

DO 20 M = MMAXH 1,MMAX 
MM = M - MMAXH 
0( L t M) = -D(L,MM) 

CONTINUE 
CONTINUE 
RF T UR N 
END 


SUBROUTINE CR2NI 

ROUTINE to SOLVE PY CYCLIC REDUCTION A BLHCK-TRI DIAGONAL 
MATRIX EQUATION WITH A NFIJMANN-LTKE CONDITION IN 2 DIMENSIONS. 

DIMENSION V ( 129 , I 30 ) , H( 7, 128 ) , HUM ( 1?9) 

COMMON NEX, NX,NEY f NY,RETSG,A, A1 t ALFA, ALFA1,K, V/DC/D 

JM = 1 + NX 
KM = 1 + NY 
LM AX = NEY 

FORWARD RECURSION — FIRST LFVEL IS L=l. 

KH « KM 
EPS = l. 

DO 1 K a 3 , KH , 2 
DO 2 J a ?,NX 

V ( J , K ) a 7 • *V ( J , K > 

CONTINUE 

ALFA 1 a A 1 
ALFA = A 
CALL GE2N1 

IF (K.EO.KM) EPS = 0. 

DO 3 J = 2, NX 

V ( J , K) a V ( J , K) + V ( J , K— 1 ) + EPS*V(J,K+1) 

CONTINUE 

CONTINUE 
DO 4 L = 2,LMAX 
NH1 a 2**(L-2) 

NH2 a 2*NH1 
NH3 = 3*NH 1 
NH4 = 4*NHl 
KL = NH4 + 1 
KH = KM 
I FPS = 1 

DO 5 K = KL, KH f NH4 

IE (K.EO.KM* IFPS = 0 
FPS a FLOATC IEPS) 

Kll a K - NH1 
K 1 2 a K + NH1*IEP$ 
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oo in 't 


K21 = K - NH2 
K 22 = K ♦ NH2+TFPS 
K3 1 = K — NH3 
K32 = K + NH3*IEPS 
DO 6 J = 2 » NX 
DUMIJ) * V ( J* K) 

VIJ,K) = 2.*VCJ,K)-V(J,K11)+Vf JfK21)-VCJ»K3l) 

1 +EPS*I-VI J,K12)+VI J,K22)-Vl J,K32) ) 

6 CONTINUE 

DO 7 M = 1 » NH 2 

ALFA I = A1 + DIL-lfM) 

ALFA = A + CIL-ltM) 

CALL GE2NI 

7 CONTINUF 
DO 8 J = 2»NX 

VI J * K } * VIJfK) + DUMIJ) - VlJtKlll ♦ VlJtK21) 
K 12) + VI J f K 22) ) 

CONTINUE 
CDN T INUE 
CONTINUE 

C++-*- BACKWARD RECURSION — DEFINE NEW INDEXt LB = LMAX-L+1 
L = LMAX 

K = KM 

NH2 = 2**1 L-l) 

NH4 = 2+NH2 

K21 = K-NW2 

DO P J = 2, NX 

DUMIJ) = V I J * K ) 

9 CON T INUE 

Dn 10 M = 1 T NH4 

ALFA1 = A1 + DlL,M) 

ALFA = A + DC L » M ) 

CALL GF2N1 

10 CONTINUF 

DO 11 J = 2, NX 

V(J»K) = 2* V I J t K) + DUMIJ) - VIJ.K21) 

11 CONTT NUE 

DO 12 LB = 2, LMAX 
L = LMAX - LB + 1 
MH2 = 2**1 L-l ) 

NH4 = 2*NH2 
KL = NH4 + 1 
KH = KM - NH4 
KI = 2*NH4 
DO 13 K = KL t KH t KI 
K21 * K — NH2 
K22 = K + NH2 
K 41 = K - NH4 
K47 = K + NH4 
DO 14 J = 2, NX 
DUMIJ) = VIJfK) 

VIJfK) = VIJfK) + V IJ,K41) + VIJ.K42) 

14 CON T !NUE 

DO 15 H * 1 f NH4 

ALFA1 = A1 + DILfM) 

ALFA * A + DIL,M) 


EPS*I-VIJ, 


NEY-L+1. 



nonofl 


15 

16 
13 
12 

13 

17 

C 

1 

c 

2 


CALL GE2N1 
CONTINUF 
on 16 J = 2 1 NX 

V(J»K) * V(J,K) + .5*(DUM(J) - V(J,K21) - V<J,K?2>) 
CONTINUE 
CONTINUE 
CONTINUE 

DO 17 K = 2tNV f 2 
DO 18 J = 2, NX 

V ( J * K) = V(JfK) + V(JtK- 1I + V(J,K+i) 

CONTINUE 
ALFA1 * Al 
ALFA = A 
CALL GF?N1 
CONTINUE 
RETURN 
END 


SUBROUTINE GE2N1 

ROUTINE TD SOLVE BY GAUSSIAN ELIMINATION the TRIDI AG. FQ • , 
TC-BETSQ,ALFA,-BETSQ)*V = F (WITH THE EI^ST OIAG. ELEMENT D E PL AC FD 
BY ALFA1), WHFPE V AND F ARE 2 -D I MENS T ON AL * 

IN THE PROGRAM, V AND F ARE EQUIVALENT, 

DIMENSION Vf 129,130) ,S(12<?) 

COMMON NEX, NX,MEY t NY, RFTSQ, A, Al ,ALFA ,ALFA1,K , V 

AL FA 1 T = 1* /ALFA 1 
BETSQI = l./BETSG 
ALFB = ALFA/BETSQ 
M2) = -AL F A II F T SQ 
V ( 2 , K 5 = V( 2 ,K ) *A LF Al T 
DO L J = 3, NX 

S( J ) = -l./CALFR + S(J-ll) 

V(JtK) = -S(J>*(BETSQI*V(J,K ) + V(J-1,KH 
CONTINUE 

FOR BACKWARD SWFEP, DEFINE NEW INDEX, JR = NX+2-J. 

NX P2 = NX + 2 
DO 2 JB = 2, NX 
J = NXP2 - J B 

V ( J , K) = V ( J , K ) - S(J)*V( J+1,K) 

CON T T N UE 

RETURN 

END 
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APPENDIX B 


FORTRAN LISTING OF EXAMPLE PROGRAM AND SUBROUTINES FOR FIRST 
REVISED CAUCHY-RIEMANN SOLVER, VERSION B 



o o o o o n 


C MAIN PROGRAM, (JBIC1B76 FOR GDC 7600 

C+*+ SOLUTION OF CAUCHY— RIEMANN FQS . FOR TNCOMDP. FLOW OVEP THIN BTCONV 
AIRFOIL USING VFRSION B OF C AUCHY-P I EMANN SOLVER. 

AIRFOIL CHORD IS FROM X=-.5 TO +.5. 

THE METHOD AND THE PROGRAM ARE DESCPIBFD IN NASA tn D-7934 
BY E. 0. MARTIN AMD H. LOMAX. 

DIMENSION U( 129, 129), VC 129,130) ,VT( 12 9) 

COMMON NFX,NX,NFY,NY,BFTSO,A, A1,ALFA, ALFA1,K,V 
C 

XJU(JU) = XL + HX*FL0AT(JU-1) 

YKUCKU) = YL + HY *f FLOAT (KU 1—1.9) 

XJVtJV) = XL + HX*(FL0AT(JV)-1.5) 

YKVCKV) = YL + HY*FL0ATCKV-1) 

RHOC X, Y) = SQRTC CCX+.5)**2 + Y*Y )/ ( ( X-. 5 ) **? ♦ Y*Y») 

ANPHI(X,Y) = A T AN?( Y ,X-.5 ) - AT AN2 C Y, X+. 5 ) 

UEX C X, Y ) = FCTR*C1. - X*ALOG( RHD C X, Y > ) - Y*ANP HI C X, Y ) ) 

VEXCX.Y) = FCTP* CY*ALOG CPHO(X,Y) )- X*ANPHI ( X, Y ) ) 

C 

FC.TP = 4. /3. 14159265 

1 F0RMAT(4I5,4F10.0,2T5) 

2 R FA DC 9, 1 ) NFX,NX,NFY, NY » XI , XU, YL , YU ,NBC , LA SC AS 
CALL SECONDCTTU 

C THE INPUT VALUES OF XL, XIJ, AND YU ARE NOMINAL VALUES, USFD TO 

C COMPUTE the EXACT VALUES. 

HY = C YU - YL) /FLOAT (NY) 

YU = YL <• HY* ( NY - .5) 

HX = tXU - XL)/ c LOAT(NX) 

JLE = 1.005 + C-.5 - XD/HX 
XL = -.5 - HX*C JLF - 1 ) 

XU = XL + HX* CNX - .5) 

JM = 1 *■ NX 
KM = 1 + NY 
NYP2 = NY + 2 
BETA = HY/HX 
PFTSO = BP?A**2 

C CALCULATE V=VTCJ> AT K=NY . 

Y = YK V( NY ) 

DO 101 J = 2, NX 
X = XJVCJ) 

IF CNBC.fo.2) GO TO 102 
VTC J) = V FX ( X , Y ) 

GO TO 101 

io? v t c j ) = o. 

101 CONTINUE 

C ++* TN I T I AL I ?E INTERIOR VALUES OF MASS SOURCFS AND VORTTC.ITY IF IS 
C FQUIV. TP V AND G IS EOUTV. TO U) — 

DO 3 K = 2, NY 
DO 4 J = 2, NX 
V C J * K ) = VTCJ) 

U(J,K) = 0. 

4 CONTINUE 

3 CONTINUE 

C++++- LOAD BOUNDARY VALUES OF U AND V — 

C IF NBC IS 1, LOAD FXACT IJ AND Vt IF NBF IS •>, LnAD U=0 AND V=0 FAR 

C FROM ATRFniL. 

K = KM 
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DO 5 J = 2» NX 
V(J,1> = 0. 

X = XJ VC J ) 

I* (X*X.LE.. 25001) V(J f n= -4.*X 

IF C ABS(X*X-.25 ).LT.1.E-51 VCJtl) = .5*VfJ f l) 

IF (NRC.EQ.l) U(JtK) * UEXfXJUf J),YU) 

IF (NBC.EQ.2) U ( J »K ) = 0* 

5 CONTINUE 

DO 6 K = ? f NY 

IF (NBC.FCU2) GO TO 33 
UCltK) = UFX CXLfYKU(K) ) 

Vf JM , K ) » VFXC XUf YKVC K) I 
GO to 6 

33 U( 1 ? K) = 0* 

V ( JM » K ) = 0. 

6 CONTINUE 

C+-M- MODIFY FRINGES OF IN T ERI0R to INCLUDE BOUNDARY VALUES — 

DO 7 J = 2» NX 

V( J » 2 ) = V f J 1 2 ) + V ( J f 1 ) 

IJ ( J ? NY ) = U( J » NY I - U C J » K M I 

7 CONTINUE 

no 9 K « 2 9 NY 

V ( 2 f K) = VC 2 ? K ) + BFTA^U(l,K) 

U C NX » K ) = UC NX t K ) + BETA*VCJM t K) 

8 CONTINUE 

C++* ZERO APPROPRIATE BOUNDARY VALUES — 

K « KM 

DO o J = l » JM 
UCJtK) = 0. 

Vf J ? K) = 0. 

VtJfl) = 0. 

9 • CONTINUE 

DO 10 K = 1 9 KM 
V f JM 9 K ) = 0. 

U f 1 9 K ) = 0. 

V(1 9 K) = 0 • 

10 CONTINUE 

C+++ NEXT DETERMINE VECTOR FO. 

DO 11 K = 2 9 NY 
DO 12 J = 2, NX 

V(J 9 K) = V(JtK) - V(J,K+1) + BETA*fU(J,K) - UfJ-l*K)> 
12 CONTINUE 

11 CONTINUE 

CALL CALCD(NEY) 

A ^ 2 • *f !• + BET$Q) 

A 1 = 2 • + BETSO 
TT2 = TT1 

T1 = SECOND(TTl) - tt? 

CALL CR2DUV 
TT2 = TT1 

T? = SECOND ( TT1 I - TT2 

C + ++ NOW VIS DETERMINED FOR J = 2 TO NX AND K = 2 TO NY. 

C NEXT DETERMINE U FOP THE SAME J AND K . 

DO 15 KB = 2 9 NY 
K = NYP2 - KB 
DO 16 J = 2 9 NX 

UCJ f KI = UCJ 9 K+II - U(J,K) ♦ BETA*(V(J 9 KI - V f J + 1 , K) ) 
16 CONTINUE 

15 CONTINUE 
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C ++4* 

c 

c 


17 


34 

IB 

19 


20 


21 

22 

23 


26 

27 


TT 2 = TT1 

T3 = SFCONDfTTl ) - TT 2 

RE LOAD BOUNDARY VALUES OF U AND V* 

IF NBC IS 1» LOAD FXACT U AND V; IF N9C IS 2t LOAD U=0 AND V=0 FAR 
FROM AIRFOIL. 

K » KM 

DO 17 J = 2, NX 

V(Jtl) * 0. 

X = XJV(J> 

IF { X*X *L F.. 25001 ) VCJtl ) = -4.*X 
IF (ABS(X*X-.25).LT.l.E-5> V(Jtl) * .5*V(J,l) 

IF CNBC.FO.l) UlJfKl = UFX( XJUC J) t YU) 

IF (NBC. EO. 2) U(J,K) = 0. 

CONTTNUF 

DO 18 K = 2, NY 

IF (NBC.FQ.2) GO TO 34 
UfltK) = UFX( XL » YKIM K1 ) 

V ( J M » K ) = VEX(XU ? YKV(K) ) 

GO TO 18 
U(ltK) = 0. 

V(JM,K) = 0. 

CONTINUE 

FORMAT ( 1 HI *$OLUT I ON OF CAUCHY-R TFMANN EQS. BY CYCLIC PFDUCTION, 
l VERS I ON B.*/ 

1* INCOMPRESSIBLE FLOW OVER T H I N RICONVFX AIRFOIL. CHORD IS FROM X= 
2-. 5 TO +.5.*//* N80=*,T3 f *. (IF NBC IS 1, EXACT BC • S ARE IMPOSED 
3FAP FROM AIRFOTL; IF NBC IS 2 t U AND V ARE 7FR0 ON OUTER BOUNDARY. 
4*// 

5* NX=*,I3t* f NY=* f I3t* f XL=*,F10.6»*, X!J=* t F 10 . 6 , * f YL=*,F10.6 f 
6* t YU=*tF10.6 f *t HX=* *E15.7,* f HY=* , Ei 5 . 7 1 * . */ /9X • 1 HJ , 3X , 1 HK , 6X , 
73HXJU, 7Xt3HYKU t 8X , 1HU, 8X , 3HUE X , 8X t 4HER RU f 1 1 X , 3H X JV , 7X , 3HYKV, 8X t 1 HV 
8,8X f 3HVEXt8Xt4HERRV) 

WRITE (6, 191 NBC f NXtNYtXLf XI),YLf YU,HXtHY 

FORMAT (6X, 2 1 4 , F 1 1 .6 , 3F 10 . 6, FI 5. 7* F 1 1. 6* 3F 1 0. 6 , E 15 .7 ) 

KI = NY/16 

IF (KI .FO.O) KI = 1 

IF ( NX -2**NFX • EQ.O ) GO TO 21 

JI = NX/20 

GO TO 22 

JI = NX/16 

IF (JI.EO.O) JT = l 

FORMAT (IX) 

DO 24 K = 1 tKM t KT 
WRITE(6t?3) 

YA = YKU(K) 

YB = YKV(K) 

DO 25 J = ltJ M tJI 
XA = XJU(J) 

XB = XJV(J) 

IF ( J.FQ.JM.OR.K.FO.l) GO TO 26 
UEX A = UEX(XAtYA) 

UA = U(J,K) 

GO T 0 27 
UFX A = 0. 

UA = 0. 

IF (J.EO.U GO TO 28 
IF (K.NF.l) GO TO 35 

TF (ABS(XP*XB-.25).GP.l.E-5) GO TO 35 
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VEXB = 0. 

GO TO 36 

3 5 VEXB = VEX ( XB ,YB ) 

36 VP = V ! J * K ) 

GO TC 29 

28 VEXB = 0. 

VB = 0. 

29 ERR U = UA - UEXA 

FRRV = VB - VEXB 

WR IT E ( 6 * 20 ) J * K* XA * YA *UA* UEX A* EPPU * XB * YB* VP* V EXP t ERR V 

2 5 CONTINUE 

24 CONTINUE 

30 FORMAT!//* T1^**E12.4** SEC,* T2=**E12.4t* SEC., T^=**E12.4** SEC. 
1*) 

WR I TE (6*30) T 1 *T 2 * T 3 

38 FORMAT!///* VELOCITIES AT Y = 0.*) 

39 FORMAT ( //26X* 1HJ * 6X, 3HXJU* IOX, 1HU,I3X ,3H UFX , 1 1 X * 3HX JV * 10X * 1HV * 13 X* 
13HVPX//) 

40 FORMAL (?3X, I4,F11.6 f 2^15.7,FIl.6*2EI5.7» 

W R I T E ( 6 * 3R ) 

WRITE! 6*39) 

00 4? J = ? * NX 
XA * XJU(J) 

XP * XJV(J) 
tIFXA - 0. 

VFXR = 0. 

IF ( ABM XA*XA-. 25) . GE.l . F-5) UFXA = UEX!XA*0.) 

IF (XR*XB.LF..2 500 1 ) VEXB * -4.*XB 
VB * V! J 1 1 > 

UA * 0. 1 2 5* f 9, *U ! J *2 ) “U ! J* 3 ) +3 •*BE’ r A* ( V ( J * 1 )-V( J+l , I ) > ) 

WRITE! 6 f 40) J, XA* UA* HEX A * XR * VR * V p XB 
42 CONTI NUF 

WR I^E! 6*45) 

DEL TX - .1 *HX 
X * 5 - I • 1*HX 

N * 0 

46 N = N + 1 

IF ! N.E0.2) X - • 5 — l • l*HX 
WR I TF ( 6* 23) 

00 47 I = 1*21 
X = X + DEL T X 
1JFXI = 0. 

VFX1 = 0. 

IF ( ABM X*X-.25 ) .GF.l.E-5) UFXl - UEX(X,0.) 

IF <X*X.L E. .25001) VFX1 * -4.*X 
WPITF(6*48) X * UEX1 * VEX I 

47 CONTINUE 

IF CN.EO.l) GO TO 46 

45 FORMAT!//* EXACT VELOCITIES AT EOGFS . */ /29X * 1HX * 1 IX * 3HU C X, l?X, 3HVF 
IX) 

48 FORMAT (2 3X * FI 0 .6 * 2E1 5 . 7 ) 

IF (LASCAS.EO.O ) GO 2 
STOP 

END 
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non 


SUBROUTINE CALCDCNEY) 

C DIMENSIONS OF D CL »M ) ARE LMAX=NEY-1 ANO MM AX= 2**LMAX. 

DIMENSION Dff>,64) 

COMMON/DC/D 
LMAX = NEY - l 
D(ltl) = SQRTC2.) 

DC 1 » 2 ) = —0(1*1 ) 

DO ?0 L = 2, LMAX 
MMAX = 2**L 
MMAXH = MMAX/2 
MMAXHI = MMAXH + 1 
DO 10 M = If MMAXH 

D(LfM) = SORT (2*+ O(L-ltM)) 

10 CONTINUE 

00 20 M = MMAXHlfMMAX 
MM = M — MMAXH 
D(LtM) = -DC L » MM I 
20 CONTINUE 

30 CONTINUE 
RETURN 
END 


SUBROUTINE CR 2DUV 

ROUTINE TO SOLVE BY CYCLIC REDUCTION A BLOCK-TR I D IA GONAL EQUATION 
FOR USE IN UBIC1B. 

DIMENSION VC129 f 130) fD{6t64 )fOUM<129) 

COMMON NFX, NX,NEY,NYfBETSQf A f A1 , ALFA, ALEAltK t V/DC/D 
C 

JM * 1 + NX 
km a l + NY 
LMAX = NEY - 1 

C+++ FORWARD RECURSION — FIRST LEVEL IS L=l. 

KH = NY - 1 
DO 1 K « 3 » KHf 2 
DO 2 J = 2 t NX 

VCJtK) = 2 **V ( J » K > 

2 CONTINUE 

A LFA 1 = A 1 
ALFA * A 
CALL GE2N1 
DO 3 J = 2, NX 

V(JfK) = V(JfK) + V (JfK-1) + VCJfK+1) 

3 CONTINUE 
1 CONTINUE 

DO 4 L = 2 T L M AX 
NHi = 2**(L-2) 

NH2 = 2*NH1 
NH3 = 3*NH1 
MH4 = 4* NHI 
KL = NH4 + 1 
KH = KM - NH4 
DO 5 K = KL t KH t NH4 
K 1 1 = K - NHI 
K 1 2 " K f NHI 


56 



oo in ^ 


I 


6 

7 

C 

0 + 4* 


14 

15 

16 
13 
12 

18 

17 


K 2 1 = K - NH2 
K22 = K + NH2 
K31 = K — NH3 
K32 = K + NH3 
DO 6 J = 2* NX 
OUM(J) = V ( J, K) 

VCJ,K) = 2.*VCJ,K)-V( J,K11)+V( J,K?1)-VC J,K31) 

1 -VC J,K12)4-V( J,K22)-V(J,K32) 

CONTINUE 

00 7 M = 1,NH2 

ALFA 1 « A1 + 0CL-1,M) 

ALFA = A 4* D(L-1,M) 

CALL GE2N1 
CONTINUF 
DO 8 J = 2tNX 

V(JtK) = V { Jt K ) + DUMCJ) - VCJ,K11) 4* V(J,K21) 

1 K 1 2) + VCJ,K2?) 

CONTINUF 

CONTINUF 

CONTINUE 

BACKWARD RECURSION — DF C TNF N^W INDEX, LR = LMAX-L+1 = N FY- 
DO 12 LR = l,LMAX 
L = NEY - LB 
NH2 « 2**CL-1) 

NH4 = 2*NH2 
KL = NH4 + 1 
KH = KM— NH4 
K I = 2*NH4 
DO 13 K = KL , KH,K I 
K21 = K - NH *> 

K22 = K + NH2 
K41 = K - NH4 
K42 = K + NH4 
DO 14 J = 2, NX 
DUM(J) = V ( J, K ) 

VC J , K 1 = VCJfKI + V ( J , K41 ) 4- V ( J , K 42 ) 

CONTINUE 

DO 15 M « l , NH4 

ALFA1 = Al + D(L,M) 

ALFA = A 4* DCL , M I 
CALL C.P2NI 
CON T INtJE 
DO 16 J = 2, NX 

V( J , K ) = V ( J , K) 4- .5*CDUMCJ) - V(J,K211 - V{ J f K?2)1 

CONT INU c 
CONTINUE 
CONTINUE 

DO 17 K = 2, NY, 2 
DO 18 J = 2, NX 

V C J , K ) = V ( J , K ) + V ( J , K— 1 ) + V(J,K+1> 

CONTINUE 
ALFA1 = Al 
ALFA = A 
CALL GF2N1 
CONTINUE 
RETURN 
END 


-VC J, 
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(■*> r> o i ~> n 


SUBROUTINE GE2N1 

+++ ROUTINE TO SOLVE BY GAUSSIAN ELIMINATION THE TRIDIAG. BQ., 

T(-RETSQ,ALEA,-BETSQ)*V = F (WITH the FIRST DIAG. ELEMENT PEPLACED 
BY ALEA1 )f WHERE V AND F ARE 2-D IMENSTONAL. 

IN THE PROGRAM, V AND F ARE EQUIVALENT. 

DIMENSION V(1?9,130),S( 129) 

TOMMON NEX , NX, NEY ,NY, BETSQt A,A1,ALFA,ALFA1,K,V 
C 

ALFA1I = 1./ALFA1 
BETSQI = 1 . /BETS Q 

ALPS = ALEA/BETJO 

S ( 2 ) = -ALFA1I*RETSQ 
V ( 2 ,K) = V ( 2 , K ) *A LFAI T 
DO 1 J = 3, NX 

S( J) = -1 ./ ( ALEE + S ( J-l ) ) 

V ( J , K) = -S( J)*(BFTSOI*V( J,K) + V(J-1,KII 

1 CONTINUE 

C FOR BACKWARD SWEEP, DEFINE NEW INDEX, JB = NX+2-J. 

N X D 2 = NX + 2 
DO ? JB = 2, NX 
J = NX P? - JB 

V ( J , K ) = V( J , K ) - S(J)*V( J+I,K) 

2 CONTINUE 
RETURN 
fn n 
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APPENDIX C 


FORTRAN LISTING OF EXAMPLE PROGRAM AND 
REVISED AND EXTENDED CAUCHY-RIEMANN 


SUBROUTINES FOR FIRST 
SOLVER, VERSION C 



o o o n o n 


C MAIN PROGRAM, UBTC1C TOR CDC 7600 

C+++ SOLUTION OF CAUCHY-R IFMANN EOS. FOR INCOMPR. FLOW OVER THIN BICONV 
AIRFOIL USING VERSION C OF CAUCHY-RI EMANN SOLVER (WITH EXTRA TERMS 
AIRFOIL CHORD IS FROM X=-.5 TO +.5. 

THE MFTHOD AND THE PROGRAM APE DFSCRI BED IN NASA TN D-7934 
BY E. 0. MARTIN AND H. LOMAX. 

DIMENSION U< 129, 129 1, V( 129,1301 ,VT(129I 
COMMON NFX, NX,NEY,NY,BFTSQ, A, A 1, ALFA, ALFA1,K, V 
COMMON/C OEP /GAMMA, DELTA ,G0, DELTA! 

C 

XJU(JU) = XL + HX+FLOATC JU-1) 

YKU(KU) = YL + HY*( FLOAT ( KU ) — 1 . 5 I 
XJVUV) = XL + HX*(FLOAT{ JVJ-1.5) 

YKV(KV) = YL + HY *FLOAT (KV— 1 ) 

RHO( X , Y) = SORT C I (X*. 5) **2 ♦ Y*Y > / ( (X-.5 ) **2 ♦ Y*Y I I 
ANPHI ( X,Y) = ATAN 2( Y* X- • 5 ) - A TAN2( Y, X+. 5) 

UFX( X , Y I = FCTR*(1. - X*ALOG( RHO (X , Y )) - Y*ANPHI ( X, Y) ) 

VE X( X, Y) = FCTR*(Y*ALOG (PHOt X,Y) )— X* ANPHI ( X,Y » } 

C 

FCTR « 4. /3. 14159265 

1 F0RMAT(4I5, 4F 10.0,1 51 

2 READ(5 , 1 ) NEX, NX ,NEY ,NY »XL,XU»YL, YU, N BC 

103 FORMAT ( 2E1 0. 0 , 15 ) 

104 PF AD (5, 103) ALPHA 1, ALPHA2,LASC AS 
CALL S ECONO ITT 1 ) 

C THE INPUT VALUES OF XL, XU, AND YU APE NOMINAL VALUES, USED TO 

C COMPU T F THE EXACT VALUES. 

HY = (YU - YL ) /FLOAT ( NY ) 

YU = YL + HY* ( NY - .5) 

HX = (XU - XL J/FLOAT ( NX ) 

JLE = 1.005 + (-.5 - X L ) / HX 
XL = -.5 - HX*( JLE - 1 » 

XU = XL + HX*(NX - .5) 

JM a 1 * NX 

KM = 1 + NY 

NYP? = NY + 2 

EFTA = HY/HX 

RETSO = BET A**2 

BETA1 = BET A*( 1 . - ALPHA1 ) 

BETA2 = BE T A*( 1 . + ALPHA2) 

BFTA^ =-RETA*ALPHAl 
BETA4 =-BETA*ALPHA2 
GAMMA = -BETA*BETA? 

DELTA = -BETA*BETA1 
GD = GAMMA/OELTA 
DELTAI = 1. /DELTA 
C CALCULATE V=VT(J) AT K=NY. 

Y = YKV(NY) 

DO 101 J = 2, NX 
X = XJVCJ) 

IF (NBC.EQ.2) GO TO 102 
VT ( J ) = VEX ( X , Y ) 

GO to 101 
102 FIJI « 0. 

101 CONTINUE 

C+++ INITIALIZE INTERIOR VALUES OF MASS SOURCES AND VORTICITY (F IS 
C EQUTV. TO V AND G IS EOUIV. TO U)~ 


60 


DO 3 K = 2, NY 
Y = YKU (K ) 

XI = XL 

IJE XI - UEX(XlfY) 

DO 4 J = 2, NX 
X = XJUCJ) 

U(JtK) = O, 

UEXX = UEX CX »Y) 

V(JfK) * VTCJ) + BFTA3*UEXX + BFTA4*UEXl 

Xi s x 

UEX1 = UEXX 

4 CONTINUE 

3 CONTINUE 

C++++ LOAD BOUNDARY VALUES OF U AND V — 

C IF NBC IS l» LOAD EXACT |J AND V; IF NBC IS LOAD U=0 AND V=0 FAR 

: FROM AIRFOIL. 

K = KM 

DO 5 J = 2 f NX 
VCJfl) = 0. 

X « XJVCJ) 

IF ( X*X.LE.. 25001) V(Jtl)= -4.*X 

I c C 4BS(X*X-.25).LT.l.E-5) V(J,1) = .5*VCJ f l) 

IF ( NBC. EQ. 1 ) U(JtK) * UEX ( X JIM J ) * YU ) 

IP CNBC.E0.2) UCJfK) * 0. 

5 CONTINUE 

DO 6 K = 2, NY 

IF (NBC.FO.?) GO TO 33 
U(1 tK) = UFX (XL » YKU ( K ) ) 

VCJM.K) = VEXCXUfYKV(K) ) 

GO TO 6 

33 UdfK) = 0. 

V(JM # K) = o. 

6 CONTINUE 

C + *+ MODIFY FRINGFS OF INTERIOR T 0 INCLUDE BOUNDARY VALUES — 

DO 7 J = 2 ? NX 

V( J » 2 ) = V f J ? 2 ) + V(Jf 1 I 
U(JfNY) * IJ(JfNY) - tJ ( J f KM) 

7 CONTINUE 

DO 8 K = 2fNY 

V ( 2 t K) = VC 2 1 K ) + BET A2*U ( 1 » K ) 

U(NX f K) = U(NXfK) + RETA*V(JM f K) 

8 CONTINUE 

C+++ 7ER0 APPROPRIATE BOUNDARY VALUES — 

K * KM 

DO 9 J = lfJM 
tJ f J t K I = 0. 

VfJ,K1 * 0. 

V(Jfl) * 0. 

9 CONTINUE 

DO 10 K = 1 f KM 
VC JM » K ) = 0. 

U(lfK) * 0. 

VClfK) = 0. 

10 CONTINUE 

C+++ NEXT DETERMINE VECTOR FO. 

DO 11 K = 2 1 NY 
DO 12 J = 2, NX 

VCJ,K) = VIJfK) - V ( J * K +1 ) + BFTAl^UCJfK) - BETA 2*U< J-l f K) 

12 CONTINUE 
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11 CON TT HUP 

CALL CALCD(NEY) 

A 1 = 2 . + RETA*BPT41 
A = A 1 + RETA*BETA2 
TT? = TT1 

T1 = SEC OND ( T T 1 ) - TT2 
CALL CR2UVC 
TT 2 = TT1 

T2 = SECOND (Till - TT2 

+++ MOW V IS DETERMINED FOR J = 2 TO NX AND K = 2 TO NY. 

NF XT DETERMINE U FOR THE SAME J AND K. 

DO 15 KB = 2, NY 
K = NYP2 - KB 
DO 16 J = ?, NX 

U { J * K ) = U ( J * K+ 1 ) - H ( J , K ) + BETA* C V ( J * K ) - V(J+1,K)) 

16 CONTINUE 
15 CONTINUE 

TT? = TT 1 

T3 = SECOND { TTl) - TT? 

+ ++ RELOAD BOUNDARY VALUES OR U AND V. 

IF NBC IS 1, LOAD EXACT U AND V? IF NBC IS 2, LOAD U= 0 AND V=0 BAR 
FROM AIRFOIL. 

K = KM 

DO 17 J = 2, NX 

V(J,1) = 0. 

X = XJV(J ) 

IF ( X*X.LE.. 25001 ) V(J,1 ) = -4.*X 

IF (ABS( X*X-.25).L T .l.E-5) v(J,l) = .**V(J,1) 

IF (NRC.EQ.l) U(JtK) = UEXIXJUI J), YU) 

IF (NBC. F0.2) U( J»K) = 0. 

17 CONTINUE 

DO ! 8 K = ? , NY 

IF ( NBC.EQ.2) DO TO 34 
U ( 1 » K ) = IJEX ( XL »YKU( K ) ) 

V(JM,K) = VEX (XU,YKV(K) ) 

GO T 0 13 

34 U ( 1 , K ) = 0. 

V ( JM » K ) = 0. 

IB CON T IN UF 

19 FORMAT (1H1 *URIC1C*//* SOLUTION OF CAIJCHY— R I FMANN EOS. (WITH EXTRA 
1TERMS) BY CYCLIC REDUCTION, VERSION C.*/ 

1* INCOMPRESSIBLE FLOW OVER T H I N BICONVEX AIPFOIL. CHORD IS FROM X= 
2-. 5 TO +.5.*//* NRC=*,I3,*. (IF NBC IS 1, EXACT BC’ S ARE IMPOSFD 
3FAR FROM AIRFOIL? IF NBC IS 2, U AND V APE ZFRO ON OUTER BOUNDARY. 
4*//* ALPHA1=*, E15.7,*, ALPHA ?=*,F 1 5. 7,*.*// 

5* NX=*,I3,*, NY=*»T3 »*, XL=*, F10 .6 , *, XU=*, E10. 6, *, YL=*,F10.6, 

6*, YU=*» F10.6,*, HX=*,F 15.7,*, HY=* , E 1 5. 7, * .* / /9X , 1 HJ ,3X , 1 HK ,6X , 
73HXJU, 7X,3HYKU,3X ,1HU,8X, 3 HUE X, 8X, 4HFRPU, 1 1 X, 3HX JV, 7X ,3HYKV, BX ,1 HV 
8,3X,3HVEX,8X,4HERRV) 

WRITE (6, 19) NBC.ALPHAl, AL PHA2 , N X,NY, XL , XU, YL,YU,HX,HY 

20 F0RMAT(6X,2I4,F11.6,3F10.6,F15.7,F11.6,3F10.6,F15.7) 

KI = NY/16 

IF (KI.FO.O) KI = 1 
IF (NX-2**NFX.EQ.O) GO TO 21 
JT = NX/20 
GO TO 2? 

21 J! = NX/16 

22 IF (JI.EO.O) JI = l 

23 RQRMAT ( IX ) 
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DO 24 K = 1 , KM,KI 
WR I TE ( 6, 23 ) 

YA = YKU(K) 

YB = YKVCK) 

DO 25 J = 

XA = XJU(J) 

XB = XJV(J) 

IF { J.EQ.JM.OR.K.EO.ll GO TO 26 
UE XA = UFX (X A ,Y A ) 

UA = U(J,K) 

GO TO 21 

26 UF XA = 0. 

UA = 0. 

27 IF CJ.FQ.l) GO Tn 28 
IF (K.NE.l) GO TO 35 

IF <AB$<XB*XB-.25>.GE.l.E-5) GO TO 35 
VFXB = 0, 

GO TO 36 

35 VFXB = VFXCXB,YB1 

36 VR = V ( J , K ) 

GO TO 29 

2 8 VFXB = 0. 

VB = 0. 

29 FRRU = UA - UFXA 

FRRV » VB - VFXB 

WR !TE<6» 20) J, K,XA f YA*UA f UEXA,ERPUf XBt YB, VB, VFXB, FRRV 
25 CON T T N UF 

24 CONTINUE 

30 FOR MA T { //* T1 =* , F 12 .4 t * SEC., T2=*,E12.4,* 5FC., T3= *, El 2. 4, * SEC. 
1*) 

WR I TE { 6 » 30 ) T 1 f T 2 »T3 

38 FORMAT C ///* VELOCITIES AT Y = 0.*) 

39 FORMAT C//26X, 1 H J , 6X,3HXJU, 10X, 1HU ,1 3X f 3HUE X , 1 1 X, 3HX J V f 10X, 1HV , 1 3 X, 
I 3H VEX / / ) 

40 FORMAT ( 2 3X, I4 f El K 6 , 2F 15 . 7 ,F 1 1 . 6 , ?F1 5 . 7 ) 

W R IT F ( 6 * 38 ) 

WR I TF ( 6 » 39 ) 

DO 42 J = 2»NX 
XA = XJU(J) 

XP = XJV(J) 
a 0. 

VFXB = 0. 

IF ( ABS(XA*XA-.25) .GE.l.E-5) UEXA = UEX(XA,0.) 

IF (XB*XR. LE.. 25001 ) VEXB = -4.*XR 
VB = VfJ,l) 

UA = 0. 12 5*(9.*!U J,2)-U( J , 3 > +3 . *BF^ A* ( V U f l ) -V < J+l , 1 M > 

WRITE (6, 40) J t XA f UA , UEXA , XB , VB , VEXB 
42 CONTINUE 

WR I T F C 6, 45 ) 

OFLTX = .1*HX 
X = — • 5 — 1. 1*HX 
N = 0 

46 N = N + I 

IF CN.E0.2) X a .5 - 1.1*HX 
WR I T E ( 6* 23 ) 

00 47 T = 1,21 
X = X + D EL T X 
U EX 1 = 0. 
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VE XI = 0. 

IF { ABS< X*X-.2F).GE.1.E-S) UEX1 = UEX(X,0.) 

IF ! X*X.L E. .25001 ) VFX1 = -4.*X 
WRITFf 6,48) X,UFXi» V p Xl 

47 CONTINUE 

IF (N.EQ.l ) GO TO 46 

45 ppRMAT ( //* EXACT VELOCITIES AT EDG>=S.*//29X ,1HX , 11X ,3HUEX, 12X, 3HVE 
IX) 

48 FORMATC23X, F10.6 ,2EI5.7) 

IF (LASCAS.EQ.O) GO TO 2 
STOP 

END 


SUBROUTINE CALCDINFYJ 

C DIMENSIONS OF OIL »M) A<»F LMAX=NFY-1 AND MM AX= ?**LMAX. 

DIMENSION Df 6 * 64) 

COMMON/DC/D 
LMAX = NEY - 1 
0(1,1) » S QRT( 2. ) 

D ( l , 2 ) = -Df 1,1 ) 
on 30 L = 2 , LMAX 
MM AX = 2**L 
MMAXH = MMAX/2 
MMAXH1 * MMAXH + 1 
DO 10 M = 1, MMAXH 

D(L,M) = SORT (2.+ 0(L~l,M) ) 

10 CONTINUE 

DO 20 M = mmaxhi, mmax 

MM = M - MMAXH 
D(L,M> = -D(L, MM) 

20 CONTINUF 

30 CONTINUE 
RETURN 
END 


SUBROUTINE CR2UVC. 

ROUTINE TO SOLVE RY CYCLIC R FDUCT inN A BLO^k— TRI DTA GONAL EQUATION 
FDR USE IN URIC 1 C. 

DIMENSION V(129, 130) ,0(6,64 ),DUMfl29) 

COMMON NEX,NX,NFY,NY,RETSO, A, A 1, ALFA, ALF A1 , K , V/DC/D 
C 

JM = 1 ♦ NX 
KM a 1 4 - NY 
LMAX = NEY - 1 

C++* FORWARD RECURSION — FIRST LEVEL TS L= 1 . 

KH * NY - 1 
DO 1 K = 3 , KH * 2 
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no 2 J = 2 , NX 

VI J *K ) = 2.*VIJ,K1 

2 CONTINUE 
ALFA! = A1 
ALFA = A 
CALL OE2N2 

DO 3 J = 2, NX 

V ( J » K ) = V ( J * K ) + V(JfK-l) + V(J,K+1) 

3 CONTINUE 

l CONTJNUF 

00 U = 2 1 L M AX 
NH1 = 2 ** C L-2 ) 

NH2 = 2*NH 1 
NH3 = 3*NH1 
NH4 = 4*NHl 
KL = NH4 + 1 
KH = KM - NH4 
DO 5 K = KL * KH » NH4 
K 1 1 = K - NH1 
K12 * K + NH1 
K 2 1 = K - NH2 
K2 2 = K + NH? 

K31 = K - NH3 
K32 = K + NH3 
DO 6 J = 2, NX 

DUN! J) = V( J,K) 

V( J * K ) = 2.*V(J,K)-V(J,K1H+V(J,K211-V(J,K31> 

1 -V(J,K12)+V(J t K22)-V(J,K32) 

6 CONTINUE 

DO 7 M = 1,NH2 

ALFA1 = A1 + D(L-lfM) 

ALFA = A + D( L-l » M ) 

CALL GE2N2 

7 CONTINUE 

DC 8 J = 2, NX 

V ( J * K ) = V ( J * K ) + DUM(J) - V(JtKll) + VIJ.K21) 

1 Kl?) + V ( J , K22 ) 

CONTINUE 
CONTINUF 
CONTI NIJF 
C 

C+++ BACKWARD RFCURSION — DEFINE NEW TNDEX, LR = LMAX-L+1 = NFY- 
DO 12 LR = 1 ? LMA X 
L = NFY - LR 
NH2 = 2** C L— 1 1 
NH4 = 2*NH2 
KL = NH4 + 1 
KH = KM - NH4 
Kl = 2*NH4 
DO 13 K = KL »KH»KT 
K21 a K - NH 2 
K 2 2 = K + NH2 
K41 a K - NH4 
K47 = K +■ NH4 
DO 14 J = 2, NX 
DUMI J) = Vf J,K) 

VI J* K ) = V { J t K ) + V ( J »K41 1 + V IJi K42 I 

14 CONTINUF 


-Vf J , 
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DO 15 M = l , NH4 

ALFA1 = A 1 +D1L,M) 

ALFA = A + 0 f L , M ) 

CALL GE2N? 

CONTINUE 
DO 16 J - 2 9 NX 

V ( J * K } = V(J,I<) 4- .5*(DUM(J) - V ( J » K21 1 - V(J,K22»> 

CONTINUE 
CONTTNUF 
COM T I N1JF 

DO 17 K = 2, MY ,2 
DO 18 J = 2, NX 

V(JtK) = V(J»K) + V ( J t K— 1 1 ♦ V(J,K + 1» 

CONTINUE 
ALFA1 = Al 
ALFA = A 
CALL GE2N2 
CONTINUE 
RETURN 
END 


SUBROUTINE GE2N2 
ROUTINE TO SOLVE BY GAUSSIAN 
7 J GAMMA, AL FA, DELTA)* V=F (BUT 
BY ALFA1), WHERE V AND F ARE 
IN _r HE PROGRAM, V AND F AP.E 


ELIMINATION the TRIDIAG. EQ., 

with THE FIRST DI AG. ELEMENT REPLACED 

2-DIMENSIONAL . 

EQUIVALENT. 


DIMENSION V<129,130) ,S(129) 

COMMON NEX,NX,NEY,NY,BFTSO,A,Al .ALFA, ALFA1,K,V 
CO MMON/COEF/ GAMMA ♦DELTA, GO, DELTA! 


ALFA1I = 1 o /ALFA l 
A» FD = A? F A*O c LT A I 
S?2) = ALFA1 r+DFLTA 
VI?, XI = V ( 2 ♦ K S *A LFA IT 
DO 1 j = 3, NX 

S(JS = l./ULFQ - GD*SIJ-1)) 

V f J ♦ K 1 = S{JS*{DELTAI*V(J,KJ- GD*V<J-1,K)) 
CONTINUE 

FOR BACKWARD SWEEP, DEFINE NEW TNDEX, JB - NX+2-J. 
NXP2 = NX 4- 2 
00 2 JR = 2 ,NX 
j = NX D 2 - JB 

V { J » K 1 = V ? J ♦ K ) - S(J)*V( J4-1,K1 
CONTINUE 

Rc T UPN 

FND 
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APPENDIX D 


FORTRAN LISTING OF EXAMPLE PROGRAM AND SUBROUTINES FOR SECOND 
REVISED AND EXTENDED CAUCHY-RIEMANN SOLVER, VERSION D 



ooooooooo 


MAIN PROGRAM, UBICID FOR CDC 7600 

+++ SOLUTION OF CAUCHY-RIFMANN EOS . FOR INCOMPR. FLOW OVER THIN BICONV 
AIRFOIL US TNG VERSION D OF CAUCHY-R IFMANN SOLVFR (WITH EXTRA TERMS 
CONTAINING CONSTANT COFFFS.). 

AIRFOIL CHORD IS FROM X=-.5 TO + .5. 

THE METHOD AND THE PROGRAM ARE DESCRIBED IN NASA TN D-7934 
RY E. D. MARTIN AND H. LOMAX. 

D IMF NS I ON U( 129,129) ,V( 129, 130 ) , VT (129 ) , VR ( l 29 ) ,UL ( 129 I, VR (129 » 
DIMENSION FKM ( 129), U T ( 129} 

COMMON NEX,NX,NEY,NY,BETSQ,A,A1,ALFA, ALFA1, K, V 
COMMON/COEF/GAMMA , DELTA ,GD,DELTAI 
C 

XJU(JU) = XL + HX*FLOAT( JU-1 ) 

YKUCKIJ ) = YL + HY*( FLOAT! KU) -1. 5) 

XJV(JV) = XL + HX*(<=LOAT( JVI-1.5) 

YKV(KV) = YL + HY*FLOAT ( KV-1 ) 

PHO ( X , Y ) = SOR T( ( { X+. 5 )**2 ♦ Y*Y) /( ( X-. 5)**2 ♦ Y*Y ) ) 

AN PH I ( X »Y) = ATAN2(Y,X-.5> - ATAM2 ( Y, X + . 5) 

UE X{ X , Y) = FCTR* ( 1. - X* ALOG ( PHO (X , Y )) - Y*ANPHI (X, Y ) I 
VEX(X»Y ) = FC-TR*( Y*AL OG (R HO { X » Y ) )- X*ANPHI ( X , Y) ) 

C 

FC T R = 4./3. 14159265 

1 FORMAT (415, 4E10. 0,215) 

2 READ ( 5 , 1 ) NEX,NX*NFY,NY,XL»XU»YLtYll,NBC, NUB DRY 

103 F0RMA T (2F10.0»I5) 

104 READ(5, 103 ) AL PHA 1, ALR HA 2, LA SC A S 
CALL SFCOND(TTl) 

C the INPUT VALUES OF XL AND XU APE NOMINAL VALUES, USED TO 

C COMPUTE THE EXACT VALUES. 

HX = (XU - xd/floatinxi 
JLF = 1.005 + (-.5 - XL ) /HX 
XL = -.5 - HX*( JLE - 1 ) 

XU = XL + H X* (NX - .51 
HY = (YU - YL ) /FL OAT ( NY ) 

JM = 1 + NX 

KM = 1 + NY 

NYP2 = NY + 2 

8ETA = HY/HX 

BETSO = BETA**2 

BETA1 = BET A* ( 1 • - ALPHA1 ) 

BETA2 = RFTA*(l. + ALPHA2 ) 

BET A3 =-BETA*ALPHAl 
BETA4 =-BETA*ALPHA2 
GAMMA = -BETA*BETA2 
DELTA = -BET A*BET A1 
GD = GAMMA/DELTA 
DELTAI = 1. /DELTA 
BETA1I = 1 . /BETA 1 

C++* INITIALIZE INTERIOR values OF MASS SOURCES and VORT I C. ITY (F is 
C EQUIV. TO V EXCEPT AT K = KM AND G IS FOUIV. TO U)-- 

DO 3 K = 2, NY 
Y = YKIJ(K) 

XI = XL 

UFX1 = UE X ( XI , Y) 

DO 4 J = 2, NX 
X = XJU(J) 

UE XX = UEX (X , Y) 
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0 


+ BFTAA*UEX1 


U ( J t K 1 = 

V ( J * K ) * R FT A3*UEXX 
XI = X 
UFX 1 - UEXX 
A CONTINUF 

3 CONTINUE 

Y = YKUCKM) 

XI = XL 

UE XI - UFXC X 1 *Y I 
DO 105 J = 2 , NX 

X = XJUCJ) 

UF XX = UE XI X t Y ) 

FKMCJ) a BFTA 3 *UEXX ♦ RETAA*UFX 1 

UC J » KMI a 0 . 

XI - X 
UEX 1 * UFXX 

105 CONTTNUF 

++++ CALC BOUNDARY VALUE 5 OF IJ AND V — 

IF NBC IS 1 , CALC EXACT U AND V? IF NBC IS 2 , CALC U =0 AND V =0 FAR 
FROM AIRFOIL, 

DO K J = 2 * JM 
VBCJ) * 0 , 

X - X J V( J ) 

IF C X*X.LE. . 2 * 001 ) VB ( J ) = ~A.*X 
I* CABSCX*X-. 25 ).LT.l.F- 5 ) VBCJ) = . 5 *VB<J) 

IF INBC.EQ.l) VT(J) ^ VEXCXt YU) 

I* f NBC. EQ * 2 ) V^CJ) * 0 * 

5 CONTTNUF 

IF INUBDRY.EO.O) GO to 116 

Y = YKUCKM) 

IF C NUBO*Y # GT, 1 ) Y a YU 
DO 117 J = 2 * NX 

IF CNBC.F 0 . 2 ) GO TO 118 
X = XJU(J) 

IJT(J) = UEX(X»Y ) 

GO TO 117 
118 U^CJ) = 0 . 

117 CONTINUF 
116 CONTINUE 

DO 6 K = 2 1 KM 

IF INBC.FQ. 2 ) GO TO 33 
UL( K) * UF X ( XL * YKUC K ) ) 

VRCK ) - VFX(XU»YKV(K) ) 

GO TO 6 

33 UL ( K ) *= 0 , 

VRCK) * 0 . 

6 CONTINUE 

C+-M- MODIFY FRINGE* OF INTERIOR TO INCLUDE BOUNDARY VALUES — 

DO 7 J = 2 , NX 

VC J ♦ 2 ) = VC J »2 ) + VBCJ) 

FKMCJ) = FKMCJ) - VTCJ) 

V C J * KM ) = FKMCJ) 

7 CONTINUE 

DO 8 K » 2 , NY 

V C 2 t K ) a V ( 2 » K ) + BET A 2 *UL CK ) 

UC NX *K) a IJ C NX f K) + RETA*VR!K) 

8 CONTINUE 

FKMC2) = FKMC2) ♦ BETA2*ULCKM1 
u « 2 v K M ) = FKMC2I 
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C++* ZERO APPROPRIATE BOUNDARY VALUES ON LEFT AND RIGHT-- 
DO 10 K = I »KM 
VCJM»KI = 0. 

U ( 1 » K ) = 0. 

V(lyK) = 0. 

10 CONTINUE 

C+++ NEXT DETERMINE VECTOR FO t FQU IV . TO V). 

DO 11 K = 2 1 NY 
DO 12 J = 2, NX 

V(J*K) = V(J,K) - V ( J t K+l ) + BET A1*U ( J y K ) - BFTA 2*U( J- 1 , K) 

17 CONTINUE 

11 CONTINUF 

C+ + * ZERO APPROPRIATE BOUNDARY VALUES ON TOP AND BOTTOM* 
on 9 J = 1,JM 

V ( J * KM ) = o. 

V<J,1) = 0. 

9 COM T INUF 

CALL CALCDtNEY) 

A1 = 2. +• BETA^BETAl 
A = Al + BE T A*RFT A2 
TT2 = TT1 

T l = SECONO(TTl) - T T 2 
CALL CR2UVC 
TT2 = TTl 

T? = RFCONJDt ^11 - T T 2 

C +++ NOW V IS DETERMINED FOR J » 2 TO NX AND K = 2 TO NY. 

C NFXT DFTERMINE U FOR ALL J AND K. 

IF (MtJBDRY.NE.O) GO To 119 
DO 106 J = 2 1 NX 

UCJtKM) = BET A1 1* C BET A2*UC J-l ♦ KM ) + V(J,NY) + FKM ( J ) ) 

106 CONTINUE 
GO TO 120 

Ho GO Tp { 121 y 122,123) NUBDRY 

121 DO 124 J = 2, NX 

UCJ»KM) = IJT ( J ) 

124 CONTINUE 
GO TO 120 

122 BETA 5 = .5*BETA 
DO 12 5 J = 2 , N X 

U ( J , KM ) = UT ( J ) + BFTA5*(VT(J) - VT(J+1)) 

125 C 0 N T T N UF 
GO ^0 120 

123 BF T A6 » .125*RFTA 
BF T A7 * 3.*BETA6 
V(JM,MY) « VR(NY) 

DO 126 J * 2, NX 

UC J, KM ) = U T ( J 1+BFTA7*< VT< J)-V T C J+l ) )+8FTA6* ( V(J, NY ) -V { J+l , NY ) ) 

126 CONTINUF 

VC JM,NY) = 0. 

120 CONTINUE 

DO 15 KB = 2, NY 
K = NYP 2 - KB 
DO 16 J = 2, NX 

U{ J, K) = U(J,K+1) - U ( J , K ) + BFTA*(V(J,K) - V<J+1,K)) 

16 CONTINUF 

15 CONTINUE 
TT2 = TTl 

T3 = SECPN D ( T T 1 ) - T T 2 
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C+++ RELOAD BOUNDARY VALUES On U AND V. 

IF NPC IS 1» LOAD EXACT U AND V; IF NRC K ?, LOAn U=0 AND V=0 FAR 
FROM AIRFOIL. 

DO 17 J = 2 ,NX 
VIJ, 1) = VB( JI 
V ( J * K M ) = VT ( J ) 

17 CONTINUE 

DO 18 K = 2, KM 
U C 1 * K ) = (JUKI 
V(JM,K) = VR ( K ) 

18 CONTINUE 

19 FORMAT (lHl*UBICin*//* SOLUTION OF C AUCHY-P I C MANN FOS. (WITH FXTRA 
ITERMS) BY CYCLIC REDUCTION, VERSION 0.*/ 

1* INCOMPRESSIBLE FLOW OVER THIN RICONVFX AIRFOTL. CHORD TS FROM X= 
2-. 5 TO +.5.*//* MBC=* , 18 . IF NPC IS 1, EXACT BC S ARE IMPOSED 

3FAR FROM AIRFOILS IF NBC IS 2, U AND V ARE Z=RO ON OUTER BOUNDARY. 
A*//* NUBDRY = *,I3,*. ROW OF U(J,K) JUST T N S IDE UPPER BOUND AR Y , T . E . 
4, U(J,KM), ts COMPUTED DIFFERENTLY ACCORDING TO SPECIFIED NUBOPY.* 
4/* TE NUPDRY= 0, U(J,KM) C OMPU t fd Fonw CONTINUITY FO. , (UNSTABLE 
4 CALCULATION I c BETA2/BFTA1 IS G T 1.0)* 

4/* = l, U(J,KM) IS COMPUTED FROM U SPECIFIED AT Y=YKU( KM ) 

4, DETERMINED FROM SAME RFLATION T HA T DFtpomtnES OUTER CONOS.* 

4/* = 2, U ( J , K M ) TS COMPUTFP USING 1ST-0P DFR-ACC . ONE-SIDE 

4D Y-DF R I V. IN ROTAT. PQ. AT Y = YU, WI T H II SPECIFIED THERF.* 

A/* = 3, Uf J , K M ) IS COMPtJTFD USING 2ND -OR DFR-ACC. ONE-SIDE 

4D Y-DERIV. IN ROTAT. FQ. AT Y = YU, WITH U SPFCIFIFD THERE. 

4*//* ALPHA1=*,E1*5.t,*, AL°HA2=**E15.7 ,*.*// 

5* NX=*, 13, *, NY=*,I3,*, XL = *,F10.6,*, XIJ=* , F 1 0. 6 ,* , YL=*,F10.6, 

6*, YU=*,F10.6,*, HX=* ,E15 .7,*, HY=*, El 5 . 7, * .*/ /9X , 1HJ , 3X , 1 HK , 6X, 
73HXJU, tx,3HYKU,8X , IHU, BX , 3HI1F X , 8X , AH ER R U, 1 1 X , 3HX J V , 7X ,3HYKV , 8X , 1 HV 
8,8X,3HVEX,RX,4HFPRV ) 

WRITE( 6, 19) NBC , NIJBDRY , ALPHA1 , ALPHA? , NX , NY, XL » XM , YL , YU , MX, HY 

20 F0PMAT(6X, 2IA, F11.6,3F10.6,F1B.7,F11.6,3F10.B,E1f.7) 

KT = MY/16 

IF (KT.FQ.O) K I = l 

IF ( NX-2**N EX . ED . 0 ) G n TO 21 

JI = NX/20 

C,n -rp ?? 

21 JT = NX/16 

22 IF (JI.EQ.O) JT = 1 

23 FQpMA-dX) 

DD ?4 K = 1 , KM, K I 
W p I TE ( 6 , 2 3 ) 

YA = YKU(K) 

YB = YKV(K) 

DO 25 J = 1,JM,J! 

XA = XJU(J) 

XB = XJV(J) 

IF (J.FQ. JM.OR.K.fq.i ) GO T 0 26 
UEXA = UFX ( X A ,Y A ) 

UA = U(J,K) 

GO to 27 

26 UFX A = 0. 

UA = 0. 

27 IF ( J.FQ.l ) GO to 28 

IE (K.NE.l) GO TO 35 

IF ( A B S ( X* X— . 25) .OE.l.E-5) RP Tp bb 
V FXB = 0. 

GO Tn 36 
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35 VEXR = VF X ( XB *YR ) 

3 6 VR1 = V ( J * K ) 
no TO 29 

28 VF XB = 0. 

V B 1= 0. 

29 EPRU = UA - UEXA 

FPRV = VB 1 - VF XB 

WP ITF(f f 20) JtKf XAtYAt UA f ueXA, FPPlJtXBt YBtVRlt VFXB *EPPV 
25 CONTINUE 

24 fON T TNUE 

30 FORMAT!//* T1=*,F12.4,* SEC.* T2=* f Fi2.4** SEC.* T3=**E12.4 t * SFC. 
1*) 

WRI T F(4*30) T1,T2*TB 
3R FORMAT!///* VFLOCT T JFS AT Y = 0.*) 

39 FORMAT l //2 6X* 1HJ * 6X, 3HXJU , 10X * 1H U ,1 3 X ,3HUFX , 1 1 X* 3HY JV , 10X , 1 wy f 13X, 
13HVTX//) 

40 FORMAT !23X* 14, FI 1 .6 * 2F 15.7 * FI ! . 6 * 2F 15.7) 

WPITF!6*39) 

DO 42 J - 2 * NX 
XA = XJIJ(J) 

XB = XJV!J) 

UF XA = 0. 

V P X B = 0. 

TF ( ARS! X A*X A— *25).GF.l.E-5) UFXA = UFX(XA,0.) 

TF ! XB *X B » L F • • 2 500 1 ) VFXB = -4.*XR 

vbi = v( j* n 

UA = 0.125*!9 # *U! J*2)-U! J* 3 )+? . *RET A*< V ! J * l)-V(J+l*l )) ) 

WR ITF! 6, 40 ) J f XA,tlA , HEX A* Xft f VBl f VFXB 
42 cnNTTMUF 

WR I T F ( 6 * 45 ) 

DFLTX - .1*HX 
X = -.5 - 1 • 1 * HX 
N = 0 

46 N = N + 1 

IF (N.F0.2) X = .5 - 1.1*HX 
WP ITF ! 6* ?3) 

PO 47 I = 1*21 
X = X + DEL T X 
U C X1 = 0. 

V C X1 = 0. 

IF <A8S! X*X-.?5).GF.l.F-5) UFX1 = UFX(X*0.) 

IF (X*X. LF.. 25001 ) VFX1 = -4.*X 
WRTT C !6*4R ) X*UFX1 * VFX1 
4T CONTIMUF 

I f CN.FO.l ) GO TO 46 

45 FORMAT!//* F X AC T V^LDC T TIFS AT F DCFS • *//29 X * 1HX * 1 IX , 3HU fx* 1 2X * 3HVF 
IX ) 

4 8 F0PMAT{?3X*F10,6* ?F15.7? 

IF (LARCAS. FQ.O) GO Tn ? 

STPD 

FND 
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SUBROUTINE CALCO ( NEY ) 

C DIMENSIONS OF DCLtM) ARE LMAX=N EY-1 AND MMAX= 2**LMAX. 

DIMFNSTON D(6,64) 

COMMON /Dr/D 
UMAX = NFY - I 
Ddtl) - SOR T C2.) 

DC 1 » 2 ) = -DC 1* n 
DO 30 L = 2 v LMAX 
MMAX = 2**L 
MMAXH = MMAX / 2 
MMAXH1 » MMAXH + 1 
DO 10 M a \ f MMAXH 

DCLtM) = SQRTC2.+ OCL-lfMI) 

10 CONTINUE 

DO 20 M = MMAXH If MM A X 
MM = M - MMAXH 
DCLtM) = -DC L t MM) 

20 CONTINUE 

30 CONTINUF 
RETURN 
END 


SUBROUTINE CR2UVC 

ROUTINE TO SOL V F RY CYCLIC REDUCTION A P LOCK-TRT DI AGONAL EQUATION 
FOR USE IN UBIC1C. 

DIMENSION V(129fI30) »D( 6» 64 ) f DUMC129) 

COMMON NEX,NX,NPYtNY,RETCO f A, A1 t ALF A, ALF A1 » K , V/DC/O 
C 

J M = 1 + NX 
KM = 1 + NY 
LMAX = N C Y - 1 
IF CNY.EQ.? ) 00 TO 20 

C + ++ FOPWAFn RECURSION — FTPFT LEVEL T R L = l. 

KH = NY - 1 
DO 1 K = 3,KH,? 

DO 2 J = 2, NX 

V ( J * K ) = 2.*V(J,K) 

2 CONTINUE 
ALF A 1 = A 1 
ALFA = A 
CALL GE2N2 

DO 3 J = 2* NX 

VCJtK) = V(JtK) + V ( J t K-l ) + VCJfK+1) 

3 CONTINUE 
1 CONTINUF 

IF (NY.EQ.4) GO TO 21 
DO 4 L = 2, LMAX 
NH1 = 2** CL-2 ) 

NH2 = 2*NH1 
NH3 = 3*NH 1 
NH4 = 4*NH1 
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KL = NHA + 1 


KH = KM - 

NHA 

no 5 K = KL , KH, 1 

K 11 - K 

- NH 1 

K 1 2 = K 

+ NH1 

K2 1 " K 

- NH2 

K22 = K 

+ NH 2 

K3 1 = K 

- NH3 

K32 = K 

♦ NH3 

DO 6 J : 

= 2 1 NX 


OUM( J) = V< J,K) 

VU,K> = ?.*VlJ,K)-VU,KllH-V(J,K?n-VIJ,K 31 ) 

1 -VIJ,K12I+Vl J,K22)-V(J,K3?) 

6 CONTINUE 

DO T M = 1,NH2 

AI.FAI = A1 + D(L-lfM) 

AL c A = A + T(L-!,M) 

CALL C=2N2 
T CONTINUE 

DO ft J = 2, NX 

V(J,K> = VJJ,K> + nUM(J) - V(JtKTl) + VIJ,K?1J 
1 K12) + V I J ,K 22 ) 

CONTINUE 
CONTINUE 
CONTINUE 
21 CONTINUE 

% * BACKWARD RECURSION — DEFINE NFW INDEX, LB = LMAX-L+l = N R Y- 
nn 12 L R •- l , LMAX 
1. = NEY - LR 
NH2 = 

NHA - ?*NH? 

KL = NHA + 1 
KM - KM — NHA 
K I = 2 *N H 4 
Or 1"> K = KL,KH,KT 
l<2 1 = K - NH? 

X 2 2 = K + NH? 

K41 = K - NHA 

KA’ = K 4- NHA 

OH 14 J = 2, NX 

niWJ) = VI J,K) 

VI J , K 5 = VI J , K ) + V I J , KA1 ) + V ( J , KA? ) 

1 4- CONTINUE 

no IS H = l, NHA 

ALFA 1 = A 1 + D(L,M) 

ALFA = A ♦ 0 1 L » M ) 

CALL OF ?N? 

VS CONTINUE 

nn 16 J = 2, NX 

V( J , K ) = VI J » K ) 4- . 5 * IDUM l J ) - V(J,K?l) - V(J,K2?1I 

16 C,nN T INUE 

13 CONTINUE 

1? continue 

20 CONTINUE 

no 17 K = 2, NY, 2 
00 IB J = 2, NX 

V I J , K ) * V I J , K) V I J * K -1 ) + VIJ,K+1) 

18 CONTINUE 


-VI J, 
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ALFA l = A 1 
ALFA = A 
CALL GF7N2 
17 CONTINUE 
RETURN 
END 


SUBROUTINE GE2N2 

C+++ ROUTINE T o ^TLVF BY GAUSSIAN FLI M IMA T IOM th= TPIDIAG. EG., 

TC GAMMA, ALFA»nFLTA)*V=F ( BUT WITH " r HF FIPS T DTAG. ELEMENT °E PLACFO 
BY ALFA1), WHFRE V AND F ARE 2— D T MENS TONAL • 

IN THE PROGRAM, V AND F ARE FQIJT VAL EN T . 

DT MENS! ON V ( 120 , 1 30 ) , S ( 1 ? P ) 

COMMON NFX , NX,NFY,NY,BFTRQ,A,A1 ,A LFA , A LF A I , K , V 
COMMON/COFF/G AR|ma ,OFLTA , GO , DF LT AT 

c 

ALFAII « 1 . / AL " A 1 
A LFD = ALFA*DELTAT 
S { 2 ) = AL r A 1 T *0E L T A 
V ( 2 , K ) = V( ?,K)*ALFA1I 
TO 1 J = 3 , NX 

S(J) = l./CALFO - G0*S(J-1>) 

V ( J » K ) = S ( j)*f CELT ai*v IJ,K )- GD*VCJ-l,Kn 

1 CONTINUE 

r FOR BACKWARD SWE C P, OF c INF NR W INDEX, JB = NX+7-J. 

NXP2 = NX + 2 
no 2 JB = 2 , NX 
J = N X p 2 - JB 

VI J , K I = V ( J , K ) - SU)*V(J + I»K) 

2 CONTINUE 
RETURN 
END 
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APPENDIX E 


FORTRAN LISTING OF 
REVISED AND EXTENDED 


EXAMPLE PROGRAM AND SUBROUTINES FOR SECOND 
C AU CHY-RI EMANN SOLVER (VARIABLE COEFFICIENTS), 
VERSION E 



oonoooo^o 


+++ 


c 


c 


c 

1 

2 

c 

r 


127 
C++ + 

c 


MAIN PROGRAM, UBIC1E FOR CDC 7600 

SOLUTION OF CAUC HY— RIFMANN EOS. FOR INCOMPR. FLOW OVER THIN BICONV 
AIRFOIL USING VERSION F 0*= CAUC HV-R I EMANN SOLVER IWITH EXTRA TERMS 
CONTAINING VARIABLE COEFFS.). 

AIRFOIL CHORD IS FROM X=-.5 TO +.5. 

THE METHOD AND THE PROGRAM ARE DESCRIBED IN NASA TN 0-7934 
BY E. D. MARTIN ANO H. LOMAX. 

DIMENSION U( 129,129) , VI 129, 130 ) , VT( 129 ) , VB 1 129 ) , UL ( 129) , VR 1129 ) 
DIMENSION UTI129) ,ALF< 12 9 ) , BET l 129 ) , GAM l 129) 

DIMENSION BETA1(129),RETA2< 129 I ,BETA3 ( 129) ,8ETA4l 129) 

COMMON NEX , NX , NEY ,NY , NX P2 , K , V 
COMMON /COEF/N,M, ALE, BET, GAM 

XJU(JU) = XL + H X *F L 0 AT (JtJ-1) 

YKUIKU) - YL + HY*IFL0AT(KU1-1.5) 

XJVIJV) = XL + HX*lFLOATlJV 1-1.5) 

YKVIKV) * YL + H Y*FLOAT ( KV— 1 ) 

RHOIX.Y) = SORT I ( I X+ . 5 ) **Z + Y*Y) /(( X-. 5) **2 ♦ Y*Y) ) 

AN PH I ( X , Y) = ATAN2I Y ,X-.5 ) - AT AN2I Y, X + .5) 

UF XI X, Y) = FCTR* I 1. - X*ALOGlRHOlX,Y) ) — Y* ANPHI f X, Y ) ) 

VEXCX,Y) = FCTR*(Y*ALOG(PHO(X,Y))- X*ANPHI t X, Y) ) 

FKJU) = XJUIJU) - .5 
F2 CJU) = XJUIJU) + .5 

FC TR = 4. /3. 14159265 
FORMAT 1415, 4E10.0, IIS) 

REAOI 5,1 ) NEX, NX , NEY , NY, X L, XU, YL, YU, NBC, NUB DRY, LA SC AS 
CALL SECOND I T T 1 ) 

THE INPUT VALUES OF XL AND XU ARE NOMINAL VALUES, USED TO 
COMPUTE THE FXACT VALUES. 

HX = IXU - XL ) /FLOATINX) 

JLF = 1.005 + I-.5 - XD/HX 
XL = -.5 - HX*I JLE - 1) 

XU = XL + H X * I N X - .5) 

HY = I YU - YD/FLOATINY) 

JM = 1 + NX 
KM = 1 + NY 
NXP2 = 2 + NX 
NYP2 = 2 + NY 
BETA = HY/HX 
BETSO = BET A**2 
ALP HA 2 = -1.0 
DO 127 J = 2, NX 
ALPHA1 = F 1 I J) 

BETA4IJ) = -BET A*ALPHA2 
BETA3IJ) » —BET A*ALPHA1 
BET A 21 J ) = BETA - BETA4IJ) 

BETA1IJ) = BETA + BETA3IJ) 

ALF(J) = BETA+BET A2 I J ) 

GAMIJ) = BETA*BETAl( J) 

BET I J ) = 2. + ALFIJ) + GAMIJ) 

ALPHA2 = F2(J) 

CONTINUE 

INITIALIZE INTERIOR VALUES OF MASS SOURCES AND VORTICITY (F IS 
EOUIV. TO V EXCEPT AT K = KM AND G IS EOUIV. TO U) — 
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4 

3 

C++ + + 

c 

c 


5 


118 

117 

116 


33 

6 

C +++ 


7 


8 

C + + + 


10 

C + + + 


DO 3 K = 2, KM 
Y = YKUCK) 

XI = XL 

UEX1 « UEXCXUY) 

DO 4 J = 2* NX 

X = XJUCJ) 

UEXX = UEX(X,Y) 

UCJ,K) a 0* 

V ( J » K ) = BFTA3C J)*UEXX + BFTA4C J )*UEX1 

XI « X 
UEX1 * UEXX 

CONTINUE 

CONTINUE 

CALC ROUNDARY VALUES OF U AND V — 

IF NBC IS 1, CALC EXACT U AND V; IF NBC IS 2, CAlC U=0 AND V 
FROM AIRFOIL. 

DO 5 J = 2,JM 
VBCJ > = 0. 

X * X JV ( J ) 

IF C X*X.LE.. 25001) VB(J)= -4.*X 

IF ( ABS(X*X-.25).LT.i.r-5) VBCJ) = .5*V8(J) 

IF CNBC. EG. 1) VTCJ) = VFX(X t YU) 

IF CNBC.F0.2) VT( J) = 0. 

continuf 

Y = YKUlKM) 

IF CNUBDRY.GT.il Y = YU 
DO 117 J = 2, NX 

TF (NBC.E0.2) GO TO 118 

x « xjum 

UTfJ) = UFX (X t Y ) 

GO TO 117 
UTC J ) = 0. 

CONTINUF 
CONTINUF 
DO 6 K = 2 t KM 

TF (NRC.FQ.2) GO x 0 33 
ULCK ) = UEXCXL, YKU(K) ) 

VPCK) a V EX C XU t YK V ( K ) ) 

GO TO 6 
ULCK) = 0. 

VRCK) = 0. 
continue 

MODIFY FRINGES OF INTER TOR TP INCLUDE BOUNDARY VALUES — 

DO 7 J = 2, NX 

V( J * 2 ) = V ( J » 2 ) + VBCJ) 

V f J » KM ) = V C J » K M) - VTCJ) 

CONTINUE 

DO 8 K = 2»KM 

V(?,K) = V C 2 ?K ) + RFTA2C2)*IJLCK ) 

UCNX*K) = UC NX » K) + BETA*VRCK) 

CONTINUE 

ZFPO APPROPRIATE BOUNDARY VALUFS ON L F p T AND 0 TGHT — 

DO 10 K = 1,KM 

V C J M f K ) = 0. 

UClfK) = 0. 

VC 1 » K ) = 0. 

CONTINUE 

NEXT DETERMINE VECTOR FO C FOU IV . TO V). 


0 FAR 
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DO 11 K = 2* NY 
DO 12 J = ^,NX 

V(J ,K) = Vf J» K) - V ( J * K + l ) + BFTAl{J>*U(J,K)-BFT A 2(J)*U(J-ltK) 

12 CONTINUE 

11 CONTINUE 

C+++ ZEPD A PPRO P RT A T E BOUNDARY VALUES ON TOP AND BOTTOM. 

DO 9 J = l*JM 
V(J*KM) « 0. 

V(Jtl) = 0. 

9 CONTINUE 

CALL CALCD(NFY) 

TT2 = TTi 

T 1 = SECOND(TTl) - TT2 
CALL CR2UVE 
TT? = TTI 

T 2 = S ECON D < TT 1 ) - t T ? 

C+++ NOW V IS DETERMINED FOP J = 2 T D NX AND K = ? TO NY. 

C NFX T DETERMINE U FOR ALL J AND K. 

119 GO TO (121* 122*123) NUBDRY 

121 DO 124 J = ? * NX 

U(J,KM) = UT(J) 

124 CONTINUE 
GO ""O 120 

12? BFTA5 = .5 *B C T A 
DO 125 J = 2 * NX 

U(J*KM| = irm + BE TA 5* ( VT { J ) - VT(J+m 

125 CONTINUE 
GO T n 120 

12? RE T A6 = .12 5#BET A 
BETA7 = 3.*BFTA6 
V ( J M , M Y ) = VP(NY) 

DO 126 J = ? * NX 

Ut J,KM) = UT(J)+BFTAT*(VT(J)-VT( J + l ) ) + RFT A6* ( V ( J * NY )-V ( J + 1 , N Y ) ) 

126 CONTINUE 

V ( J M * NY ) = 0. 

120 CONTINUE 

DO 15 KR = 2* NY 
K = NY P? - KB 
DO 16 J = 2, NX 

U ( J * K I = U(J»K+1I - !J(J,K) + RETA*(VU,K) - V(J+1*K>) 

16 CONTINUE 

15 CONTINUE 

TT 2 = TTI 

T3 = SECOND ( T T1 ) - TT? 

C+ + * RF LOAD BOUND ARY VALUF C OF II AND V. 

IF NBC IS 1* LOAD c XACT !J AND V; IE NBC IS 2* LOAD U=0 AND V = 0 EAR 
FROM AIRFOIL. 

DO IT J = 2 * NX 
V(Jtl) * VB ( J ) 

V(J ? XM) - VT ( J ) 

17 CONTINUE 

PO 18 K = 2,KM 
U( l *K) = UL ( K) 

V(JM*K> = VR(K) 

18 CONTINUE 

19 FORMAT (1H1*UBIC1F*//* SOLUTION OF CAUCHY-R I EMAMN EOS. (WITH FX T& A 
1 TERMS AND VARIABLE CO c FFS.) BY CYCLIC REDUCTION, VERSION F . * / 

1* INCOMPRESSIBLE FLOW OVER THIN BICONVEX AIREQTL. CHORD IS FROM X= 


79 



2—. 5 Tn + . 5 . *//* NBC=*,I3,*. IP NBC IS It EXACT BC' S ARE IMPOSFD 
3FAR FROM AIRFOIL; IF NBC IS 2t U AND V ARE ZERO ON OUTER BOUNDARY. 
4*//* NUBDRY=*,!3,*. ROW OF IJIJ,K| JUST INSIDE UPPER BOUNDARY, I. E 
4, IK J, KM), IS COMPUTED DIFFERENTLY ACCORDING. TO SPECIFIED NUBDRY. 
4/* TF NU BDR Y= I, IJ(J,KM) IS COMPUTED FROM U SPECIFIED AT Y=YKU<KM) 

4, determined FROM same relation THAT determines OUTER CONDS.* 

4/* = 2, U(JtKM) IS COMPUTED US TNG IST-ORDER-ACC . ONE-SIDE 

4D Y-DERTV. TN ROT AT . ED. AT Y = YU, WITH U SPECIFIED THERE.* 

4/* = 3, U( J» KM ) IS COMPUTED USING 2ND— ORDE R- ACC. ONE-SIDE 

40 Y-DERIV. IN ROTAT. EQ. AT Y = YU, WTTH U SPECIFIED THERE.*// 

5* NX= *» I 3, * , NY=* ,13,*, XL=* , FI 0.6,*, XU=* , F10 .6 ,*, YL=*,F10.6, 

6*, YIJ=*,F10 .6, *, HX=*,E15.7,*, HY=*,E15.7»*.*//9X,1HJ ,3X,1HK,6X, 
73HXJU,7X,3HYKU,8X,1HU,BX,3HUEX,8X,4HERRU,11X,3HXJV,7X,3HYKV,8X, 1HV 
8, 8X, 3HVEX, 8X, 4HERRV ) 

WRITE (6, 19) NBC , NUBDRY , NX , NY, XL, XU, YL , YU, HX , HY 

20 FDRMAT(6X,2I4»F11.6»3F10«6»E15.7»F11.6»3F10.6,E15.7) 

KI = NY/16 

IF (KI.EO.O) KI = 1 
TF INX-2**NEX.E0.0) GO TO 21 
JI = NX/20 
GO TO 22 

21 JI = NX/16 

22 IP (JI.EQ.O) JI = 1 

23 FORMAT ( IX) 

DO 24 K = 1,KM,KI 
WR I T F (6 , 23 ) 

YA = YKUCK) 

YR = YKV ( K ) 

DP 25 J = 1 , JM, JT 
XA = XJUCJ) 

XB = XJVIJ ) 

IF (J.EQ. JM.OR.K.EO.l ) GO TO 26 
UF XA = (1EX( XA ,YA ) 

UA = U( J ,K ) 

GO TO 27 

26 UF XA = 0. 

UA = 0. 

27 TF ( J.EO.l ) GO TO 28 
IF (K.NE.l) GO T 0 35 

IF <ABS(X*X-.25).GE.l.E-5) GO TO 35 
VF X° = 0. 

GO TO 36 

35 VFXB = VFX ( X R ,Y B ) 

36 VB1 = V ( J , K) 

GO to 29 

28 VFXB = 0. 

VB1 = 0. 

29 FRRU = UA - UEXA 

FRRV = VB1 - VFXB 

WRT T EC6,20) J , K , X A , YA ,UA ,UFX A , ERRU ,XB, YB, VBl , VEXB, ERPV 
25 CONTINUE 

24 CONTINUE 

30 FORMATf//* ti = *,fi2.4,* SEC., T2=*,E12.4,* SEC., T3=*,E12.4,* SEC. 
1 *) 

WRITE ( 6, 30) T1,T2,T3 

38 FORMAT!///* VELOCITIES AT V = 0. *) 

39 FORMAT ( //26X ,1HJ, 6X.3HXJU, 10X, IHU, 13X, 3HUEX » I1X, 3HX JV , 10X, IHV, 1 3X, 
13HVEX//) 

40 FORMAT 1 23X, I4,F11.6,2E1*5.7n,fii.6.2E1f.7> 
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WRITF(6,38) 

WP.ITFC 6, 3*?) 
no 42 J = 2 , NX 
XA = XJU(J) 

X8 = XJV(J) 

U EX A = 0, 

VF X8 = 0. 

IP (ABSC XA*XA-»25),GF*l*E-5) UFXA = UFX(XA,0.) 

I c (XB*XB.LE.. 25001) VEXB = -4**XB 
VP 1 = V(J,1) 

UA = 0.1?5*(<>.*U( J,2)-U( J,3)+3**RFTA*( V< J,1)-V(J+1,1 > ) ) 

WRITE 16, 40 ) J,XA,tJA,UFXA,XB, VB1, VFXB 
42 CONTINUE 

WR TT F ( 6, 45 ) 

OELTX = • 1 * HX 
X = -.5 - 1 • 1 *HX 
N = 0 

46 N = M + I 

IF ( N • C Q • 2 ) X = .5 - 1.1*HX 
WRTTP(6,?3) 

DO 47 T = 1,21 
X = X + D FL T X 
UEX 1 = 0. 

VF X 1 = 0. 

T r (ABE(X*X-.2F).GF.l.E-5) UE X l = UFX(X,0.) 

IF ( X*X*LF. .25001 ) V c Xl = -4 . *X 
WPITF(6,4R) X,UFX1, VFXl 

47 CONTINUE 

IF (N.FQ.l ) rp TP 46 

45 FORMAT!//* p XAC T VPLnCT T IFS AT F DOF c • * / /? q X , 1HX,11X ,3HUFX,12X,3HVF 
IX) 

4 8 F0PMA t (?3X,F10*6 ,2E15.7) 

IF (LASCAS.FO.O) GO t P 2 

STPo 

FND 


SURRPUTTNF CAL C0( NFY ) 

r 0 1 M r N R I PM 5 OF D ( L , M) A D F L MA X= NF Y-l AMP MMAX = ?**LMAX. 

DIMENSION 0(6,64) 

common/dc/d 

L M A X * NEY - 1 

0(1,1) = S OPT ( 2 * ) 

P( 1,2) = -0(1,1) 

DO 30 L = 2 , LM A X 
MMAX = 2**L 
MMAXH = MMAX / 2 
MMAXH1 = MMAXH + 1 
00 10 M * 1, MMAXH 

0 ( L , M) = PORT(?.+ 0 ( L-l , M ) ) 

10 CONTINUE 

DO 20 M = MMAXH1,MMAX 
MM = M - MMAXH 
0(L,M) = - D( t , MM ) 

20 CPNTTNUF 
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30 C0M T IN UE 
RETURN 
END 


SUBROUTINE CR2UVE 

ROUTINE TO SOLVE BY CYCLIC REDUCTION A BLOCK-T Q 1 0 1 A HON AL 
MATRIX EOUATION WITH VARIABLE CPE^FS. C 0R USE IN UBIC1E. 

DIMENSION V( 129,130),ALF(129) , RFT ( l?Q ) f r,AM( 129 ) ,0IJM(129> 
COMMON NEXt NX, NFYtNYt NXP2t K f V 
COMMON /CDF F /N, M f ALE, BET f GAM 

JM = 1+NX 
KM = 1+NY 
L2MAX = MFY-1 

FORWARD RECURSION — FIRST LFVEL TS L 2= I 

KH = NY-1 
DO 5 K = 3*KH,2 
DP 6 J = 2*NX 

V c J , K ) = ? . *V ( J , K ) 

CONTINUE 
N = 0 

M - l 

CALL GFIJV C 
D p 7 J = 2 * N X 

= V(J,K) + VfJ,K-l) + VCJfK+1) 

cpnttnuf 
b CONTINUE 

no 8 L2 = 2 ? L2 M A X 
N = L? - I 
MH1 = ?**(l?-2) 

MM2 = ?*NH1 
N mi = 3*MH1 
NH4 = 4 * N H 1 
KL = NW4 * 1 


KH = NY 

- NH4 

+ 1 

D n 9 K = 

: KL » KH 

1 9 NH4 

Kll * 

K-NHl 


K 1 2 = 

K+NH1 


K 2 1 - 

K-NH2 


K 2 2 = 

K + NH2 


K3 1 = 

K-NH3 


K32 = 

K+NH3 


no io 

J = ?, 

NX 


HUM I J) * V( JtK) 

V(J,K) = 2.*V( J,K)-VU,K11 )-V( J,Kl?)+V(J,K2i)+V (J, K22) 
1 -V(J,K31 )-V( J,K32) 

10 CONTINUF 

DO U M = I f NH? 

CALL GEUVE 

11 CONTINUE 
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o o xO non 


no 12 J = 2 * NX 

V{ J,K) = V( J,K) + DUMf J >-V (J »Kll)-vr J,K1?)*V< J,K211 *rVf J,K?n 
12 CONTTNUF 

9 CONTINUE 

8 CONTINUE 

BACKWARD RECURSION — DEFINE NEW INDEX? L2R = L2MAX-L 2 + i =NE Y-L? - 

DO 21 L2R = 1 , L2M AX 
12 = NEY-L2B 
N = L2 

NH 2 = 2** ( L 2- 1 ) 

NH4 = 2 * N H ? 

KL « NH4 + 1 
KH = KM-NH4 
K I = 2 *NH4 
DO 22 K = KL,KH,KT 
K 2 1 = K-NH2 
K22 = K+NH2 
K41 = K-NH4 
K42 = K+NH4 
DO 23 J = 2, NX 
OUM { J ) * VI J,K) 

V ( J, K ) = V ( J , K 41 ) + V I J , K 42 ) + V<J,K) 

23 CONTINUE 

DO 24 M = 1 , NH4 
CALL GEUVF 

24 CONTINUF 

on ?s j * 2 ♦ nx 

V { J t K ) = V ( J? K ) + .5*(DUM(J )-Vf J,K21 )-V< J,K22)) 

25 C ONTTNUF 

22 CONTINUE 

21 CONTINUE 

DO 26 K = 2 , NY t 2 
DO 27 J = 2, N X 

V ( J , K ) = VIJ,K-1) + V(JtK + l) + V ( J , K ) 

27 CONTINUE 

N = 0 

M = 1 

CALL OEUVF 

26 CONTINUF 
RETURN 
END 


SUBROUTINE GFUVE 

ROU^INF to SOLVE BY GAUSSIAN ELIMINATION THE T RT D I AG ® FQ., 
T(— ALF(j ) f BETP CJ) t-GAM(J) ) *V = F, WHF R F V AND F APF 2-D* 

DIMENS TON V{ 129, 1 30) , ALF ( 129 I iBFT< 1291 ,GA M I 129 ) , 

1 S C 129) ,016,641 

COMMON NEX , NX, NEY , NY , MXP2 ,K , V 
COMMON/DC/ 0 /CO EF/N , M , ALF f BET, GA M 
C 

S( 1) = 0* 

IF (N.EQ.O) GO TO 3 
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LL» 1 J = 2, NX 

BFTP = BE T IJ» + D<N,M> 

T = l./CBETP + ALF(J )*S{ J-l) ) 

S(J) = -GAM{J)*T 

V ( J » K ) = T*(V{J,KI + A!_F(J)*V( J-l, KM 

1 CONTINUE 
GO TO 4 

3 DO 5 J * 2, NX 

BFTP = BETCJ) 

T = 1* /( BE T P + ALF{ J)*StJ-l) ) 

S(J) = -G AM ( J ) *T 

VIJ,K) = T* ( V ( J »K ) + ALF ( J ) *V ( J-l » K ) ) 

5 CONTINUE 

C FOR BACKWARD SWEEP, DEFINE NEW TNDFX, JR 

4 DO 2 JB = 2, NX 

j = NXP2 - JB 

V ( J , K ) = V( J »K ) - S(J)*V(J + 1»K) 

2 CONTINUE 
RETURN 
END 


NX+2-J. 
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APPENDIX F 


CYCLIC REDUCTION FOR VERSIONS B TO E OF CAUCHY-RIEMANN SOLVER 


A detailed algorithm is outlined here for the recursive cyclic reduction 
of the block-tridiagonal matrix equation in the form of equations (50) and 
(59) for use in the FORTRAN programs for versions B to E of the extended 
Cauchy-Riemann solvers. The algorithm is a variant of Buneman’s (ref. 2) 
double cyclic reduction in which scalar-tridiagonal-equation solutions are 
obtained by the Thomas algorithm (Gaussian elimination). Recursive cyclic 
reduction was devised by Prof. G. Golub with collaboration of Dr. R. Hockney 
(see refs. 1 to 3). The method is based on odd/even reduction, which was 
used extensively by Hockney (ref. 1) in direct two-dimensional Poisson solvers 
and was the basis for the extension by Buneman (ref. 2) to his double cyclic 
reduction algorithm for solving Poisson T s equation. Refer to reference 5 for 
a treatise on the development of the method. 


The block-tridiagonal matrix equation to be solved represents the system 
of n matrix equations: 

-Vl + c o\ - Vi ■ »fc 0) 1, 2, . . n) (FI) 

where n = n.\ - 1, n\ = 2-k with L an integer, and in which each is 

a known m-dimensional vector. 


P 


( 0 ) 

k 



col[F l,?<’ 




Each V k is an m- dimensional vector to be determined: 


(F2) 


\ = col[y ! X 


V 7 ] 

m,k 


where, for simplicity, we have set, in equation (FI), 


(F3) 


and Cq is the m - dimensional tridiagonal matrix C defined by equa- 
tions (21c) and (21d) : 


(F4) 


C 0 Y ') 


m 


0 j- 


(F5) 


with aj, and yj defined in terms -of the parameter B and in terms of 

and whicn are specified for j = 1 to m in each version of the 

solver. 
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The cyclic-reduction algorithm depends on the matrix factorization defined 
as follows: each level of the recursion process is denoted by integer £, 
which ranges from 1 to L-l. For each £, define an integer N and quantities 
d ^ as follows: 


N = 


£, 1 


d 


a,M 


2 (A - 1, 2, ...» L - 1) 

^ ^,2--^,! « = 1 > 
^2 + dl . ,, 


d i,M d SL,M-hN 


/ £ — 2 , 3 , . • • 9 L 1 A 
\M = 1, 2, . . j N ) 

/a = 2, 3, . . . , L - 1, \ 

+ 1, . . N - 1, N/ 


(F6) 

(F7a) 

(F7b) 

(F7c) 


Then, given C Q , for each 5, there is a matrix defined by 


C 

£ 


= C? 

£.-1 


21 


(F8) 


which has N tridiagonal matrices as factors: 


(1) (2) (N ) 

C £ C £ 


(F9) 


in which the Mth factor is 


cf> = C + d. I 

£ £,M 


( £ = 1 to L - lA 
M = 1 to N / 


and 


= T (-a . , 3 ,-Y .) 

m' 3 j’ o 


3 f . = 6 . + d n 

0 3 


(FlOa) 

(FlOb) 

(FlOc) 


where -aj, 3j, and -yj are the elements of C 0 in equations (F5). (Note 
that the values of can a ^ so be obtained from a cosine function (e.g., 

eq. (21b) in ref. 25), but are calculated in the sample FORTRAN programs 
listed in appendices A to E by use of equations (F7). Note also the differ- 
ence in notation, where here has the same role as C $ in ref. 25.) 

For the first level of reduction (£ = 1), multiply each of the even 
equations (k *= 2, 4, . . . , n - 1) in equations (FI) by Cq and add to it 
the adjacent equations above and below to obtain the reduced system 
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-V, + C V, - V, L - r 7 (l) (k = 2, 4, 

k-2 i fc k + 2 k y * ’ 


. . , n - 1) 


where 


= P; 


(0) _ „(o) 

k 

, (i) 


= 1 , 2 , . 




c n r<°> +f <°) +r <°> 

o k k- 1 k+i 


. , n) 

(k = 2, 4, 


, n - 1) 


(Fll) 


(F12a) 

(F12b) 


A special device was introduced by Buneman (ref. 2) for avoiding all matrix 
multiplications (and a consequent stability problem) in the final solution 
process. The device removes the matrix Cq from the right side of equa- 
tion (F12b) and replaces it by so that the implied matrix multiplication 

need never be performed. To use this device, since we wish to remove Cq 
and insert C ls and since Cj contains C^, we write equation (F12b) as 


- 'i[‘W > ] + •£! 

- < C S - 2»[V'* <0) ] + *[«W > ] + + 


(o) 

fefi 


Now define and by denoting the first term by C^q^ ^ and the 

remainder by p^ \ Thus 


( 1 ) 


O),. „(U 


where 


■* -= c iH /+ *k 

q O) = c-ip(°> 

q k 0 P k 

( 1 ) 9n (1) . „(0) , (0) 

p k = 2q /c + v k - 1 + p k+i 


(F13a) 

(F13b) 

(F13c) 


Note that equation (F13b) is just a tridiagonal matrix equation that can be 
solved for q^ 1 ^, and then p^ 1 ^ is obtained from equation (F13c) . Also, 
substitution of equation (F13b) into (F13c) gives an equation that can be 
solved for without knowing 



For the second and higher levels of reduction (£ = 2, 3, 
we let 


(F13d) 
L - 1), 


h = (2) l ~ 2 = N / 4 


(F14) 
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and in the remaining reduced set of equations we multiply each of the even 
equations ( k = N , 227, 3 27, . . . , n\ - 27) by and add to it the adjacent 

equations above and below to obtain, for £ = 2 to L - 1: 


-V* + V* - Va- ■ 'k* (k - »• =»• » - « (F15) 


where 


r oo_ r -u-i) . r a-D . r a-D 

k - M k k-2h k+2h 

(a) (£) 

As in equations (F13) above, we wish to define q^, and so that 


(F16) 


'* l) 5 + ^ l) 


(F17a) 


and this can be done by first substituting equation (F17a) (treating as 

a function of its superscript £ and its subscript k) into each term on the 
right side of equation (F16) and then using a procedure directly analogous to 
that used to obtain equation (F13) from (F12) to eliminate C„ ^ in favor of 
C £ . We then have equation (F17a) in place of (F16) , with the results 


q 


(£) 

k 


(£-0 


+ 



(£-D 

k 


+ q 


(£- 1 ) 

k-2h 


+ q 


(£-i) 

k+2h 


P 


(£) 

k 



+ p 


(£-D 

k-2h 


+ P 


(£-0 

k+2h 


(F17b) 

(F17c) 


P K 


Furthermore, one can obtain (analogous to eq. (F13d)) an equation to solve for 

(£) 

without knowing any of the q£ . If equation (F17c) is- solved for 

(£) K (£) 

in terms of functions p£ , and the resulting q^ is then treated as 
a function of its subscript k and superscript £ to obtain expressions for 
(£- 1 ) (£- 1 ) , (£- 1 ) 
q k ’ q k-2h ’ and q k+2h 


(£) 
k 

4 


(•'.) 


to substitute into (F17b), one obtains 


(£) _ „ (£-0 
k-2h 


= P 


+ P 


(£-0 

k+2h 


+ P 


(£-0 


- P 


(£- 2 ) 

k-h 


, (A-2) 
k+h 




(£-0 


- P 


,(£- 2 ) 

k-3h 


- P 


(£- 2 ) 

k+3h 


] 


(£- 2 ) 

k-h 


,(£- 2 ) 

k+h 


+ P 


(£-0 

k-lh 


+ P 


(£- 1 ) 

k+2h 


(F17d) 


The procedure for the algorithm to solve equations (FI) , to determine all 
uses the equations derived above as follows: 

1. First, for integer L = logpfti, compute the array d^ 9 M f° r £ = 1 
to L - 1 and M - 1 to if using equations (F6) and (F7). 
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2. For £ = 1, start the forward recursion with given by equa- 
tion (F2) for k = 1 to n. Use the Thomas algorithm (Gaussian elimination) 
to solve equation (F13d) in the form 


Q Q O ) — 2d ^ ^ 

L o°k “ 2p k 


(F18a) 


successively for k = 2, 4, 


defined in equation (F5) and where 0^ ' is defined by 


n-1, where C Q is the tridiagonal matrix 

(D,, 


, (1) , 


( 0 ) 


( 0 ) 


+ »k- i + p k + 1 


After solving equation (F18a) for each unknown vector 0?: 

( 1 ) ^ 

from (F18b). 


( 1 ) 


(F18b) 

then determine 


3. For each £ = 2, 3, 


form 


. , L-ly first write equation (F17d) in the 




(i-2) 

’ k-h 


a- 2 ) 

’ k+h 


a-i) 

’k-2h 


U-i) 

'k+zh 


- P 


U"2) 

k-3h 


- P 


U-2) 

k+3h 


for each k - l V , 2/1/, . . . , n\-N , where 0 


U) 


is defined by 


(F19a) 


00 

'k 


XD 

} k 


U-i) 


(£- 2 ) 


. 01-2) 
fc+/z 


(^-D 




+ v k+2h 


(F19b) 


Equation (F19a) is of the form (cf. eq. (F9 ) ) : 


Xl)[ c (2) r (y) fl (£)1 _ Z (D 

• • • S-1 % J " k 


where 


( 1 ) 


(F20a) 


and each C 


represents the total vector on the right side of equation (F19a) 

a tridiagonal matrix defined by (F10) . Similarly, let the 

bracketed quantity on the left side of equation (F20a) be denoted by 
so that 


(2)[ (3) r W ft (£)l _ 

VlLS-1 * * * J “ 


(2) 


(F20b) 


(2) 


etc. Thus the procedure at this step is to first consider z £ J as the 
unknown vector in equation (F20a) , and use the Thomas algorithm to solve (F20a) 
( 2 ) 


for 


’k 


, then with that vector known, solve equation (F20b) for the brac- 


(£) (£) 

keted quantity, and so forth, until 0^ is known. Thus has been 

obtained by a sequence of tridiagonal solutions, and then I s obtained 
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from equation (F19b). When this is done for all k - N 9 2 N 9 . . n\-N 9 
increase Z and repeat this step 3. 

4. After step 3 has been done for Z - L- 1, all the p£ are known 

that are required for the backward recursion (back substitution). The back 
substition proceeds successively, for each l from L-l to 0. 

For Z = L- 1, L- 2, . . 1, use equations (F15) and (F17a) (with (Fll) and 

(F13a) as well) in the form 


for k = N, 3 N, 5/7, 


c & = V + y + p 

JTfe k-N k+N V k 

(Z) 

. , ni~N, where is defined by 

,<* 

l k 


V — A (O i n (O 


(F21a) 


(F21b) 


or with equations (F17c) and (F13c) , 


\ ’ *k l) - 


1 

2 [' 


a) _ u-D_ a- 

v k v k-2h p k+2h 


r] 


(F21c) 


for k = /7, 3 N 9 5N 9 . . ., Thus, one can now solve equation (F21a) , by 

use of the Thomas algorithm for a sequence of tridiagonal solutions in a 

manner analogous to that used in equations (F20) , to obtain <j>^ , with which 

V^, is then obtained from (F21c) for k = N 9 3/7 , 5/7, . . ., 

5. For Z = 0, equation (FI) can be written as 


C.v 7 = v 7 + V 7 , 4- P 7 

o k k - l k + l k 


(o) 


(F22) 


and solved using the Thomas algorithm to obtain for k = 1, 3, 5 , . . . , n. 

In a machine computation using this algorithm, the components (as needed) 

(Z) 

of each of the vectors p^ , and can all occupy the same array in 

computer memory for one value of k . Therefore, in addition to a relatively 

small array for d £ $ given by equations (F7) , only one large array, 

(tt?i+ 1) x (m+1) , is needed (Vj y) along with an m\ -dimensional dummy vector for 

5 (£) (£) 
the intermediate storage either of in equations (F18) or of <j>^ in 

equations (F21). In addition, the Thomas algorithm itself requires one dummy 

mi vector. 
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